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BY 
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BAKERY HOUSE, 40 HIGH ST., WIMBLEOON, LONDON, SWl 9 5AU 

ENGLAND 


INTRODUCTION 


1.1 Background 

The SHIP (Supersonic Hydrogen J[njection Program) computer 
program is concerned with the numerical computation of three- 
dimensional flow situations arising when hydrogen is injected 
into a supersonic airstream at an angle ranging from normal 
to parallel to the airstream main flow direction. The flow 
may be free, or confined in a duct whose walls may expand or 
contract arbitrarily but smoothly. 


In its present version, SHIP predicts the mass fractions of 
the products of combustion . of hydrogen and oxygen based on 
the assumption of chemical equilibrium. 


1. 2 Connections with Previous Work 

Under a previous Contract (NASl-13239) a computer program, 

HISS, was developed for similar situations to those considered 
here. The code has been extended and extensively modified to 
become SHIP so that the following facilities, not provided 
for in HiSS, are included: 

(i) Any of the four lateral boundaries can be (1) a wall, 

(2) a symmetry plane, (3) a free surface. 

(ii) For walls, the distance of each wall from a reference 
plane may be specified as an arbitrary, but smooth, 
function of distance along the principal flow direction. 

(iii) Ability to specify mass flux through the top and 
bottom wall is provided. 

(iv) Choice between specification of an adiabatic wall or 
a constant temperature on each wall is provided. 

1 . 3 Purpose of the Present Report 

The terms of the NASA contract call for the development and 
jbransmission of the SHIP code having the general features 
described above. 


This report marks the completion of the project. It provides 
all the necessary information concerning the mathematical 
modelling of the flows under consideration, it describes the 
numerical analysis involved in the solution of the relevant 
equations and it defines the user-orientated parts of the 
code • 

It is to serve as both a comprehensive reference to the 
mathematics and numerical procedure used in the code and as 
an operational manual for the computer program SHIP. 

1 . 4 The Physical Problem Considered 

NASA is active in research on supersonic combustion {2-4}. 

Of particular interest are methods for injecting hydrogen in 
a supersonic airstream in a manner which optimizes the design 
of supersonic combustors. 


Three arrangements of interest are shown in Figures (1), (2) 
and (3). 

Due to interaction of the mainstream with the jet, the flow 
separates ahead of the jet and reattaches behind it. Also a 
bow shock is caused by this interaction. 

There is an interaction between the jets causing flow in 
directions normal to the main-stream direction. This 
is the phenomenon which makes three-dimensional calculations 
necessary. 

Downstream of the jet the hydrogen mixes with the air and the 
region containing hydrogen widens. 

As the flow moves downstream, the mixing, chemical reaction, 
acceleration, viscous effects etc. all combine to produce a 
pressure variation. 

Thus, in addition to the aero- and thermodynamic characteristics 
of these flow arrangements, knowledge of the extent of the 
region of hydrogen, and its distribution within it, is necessary 
to assess the effectiveness and uniformity of air-hydrogen 
mixing. 


1 . 5 Capabilities and Limitations of SHIP 

SHIP is a general, flexible computer program capable of 
calculating three-dimensional, boundary- layer flows which are 
either external or internal. Features built into this 
program include: 
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FIGURE 1: DEFINITION OF GEOMETRY AND PLANES OF PROFILES FOR CASE 1 



(TYP) 


FIGURE 2: DEFINITION OF GEOMETRY AND PLANES OF PROFILES FOR CASE 9 
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NOTE: All dimensions in meters 



NOTE: Calculations are performed between 
boundaries shown in the view at the 
riqht. 




Case 7 


FIGURE 3: DEFINITION OF CASE 7 TO DEMONSTRATE THE VALIDITY OF 

THE SHIP CODE 
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• Any of the four boundai'ies cau be either a i^all, a 
symmetry plane, or a free surface. 

• For internal flows the distance of the duct walls from a 
reference plane may be specified as a function of distance 
along the principal flow direction. 

• Walls can be specified as being either adiabatic or at 
constant temperature. 

• Injection of hydrogen at angles ranging from normal to 
parallel to the main air-stream can be handled. 

• Injection through arrays of holes in either or both the 
top and bottom walls is provided for. 

• Free-boundary conditions resulting from a small-disturbance 
theory and isentropic considerations are applied to 
supersonic free boundaries. 

Thermodynamic equilibrium is supposed to prevail between the 
species Hg, 0^, HgO, 0, H and OH; and four equilibrium 
reactions are'^allowed for. 

The upper surface of the domain of integration is taken above 
the bow shock so that the edge conditions may be the same 
as the free-stream conditions. 

The lateral and longitudinal pressure gradients are uncoupled 
in the subsonic-flow in order that the equations remain 
parabolic, and forward-marching integration can be used; 
thus, the calculation "jumps" the separated region at the 
jet exit. 

The practice used for this jump is as follows: 

1. An arbitrary long forward step is taken over the 
jet exit [e.g., 10 jet diameters]. This effectively 
jumps the elliptic region near the jet exit. 

2. The source and sink terms for all of the dependent 
variables are modified at the cell immediately 
over the jet exit, so as to account for the 
required inflow rate at the jet boundary of the 
variable in question. The overall conservation 

of mass, momentum, and energy, is thus satisfied 
for the region over the jet. 

The flow is considered to be turbulent and the effective 
viscosity is calculated by way of the "k-e" two equation 
turbulence model described in {5}. Density is calculated 
via the ideal-gas law. 

The program predicts the distributions of the following 
variables : 

(i) the three velocity components; 

(ii) static pressure; 

(iii) total enthalpy; 

(iv) mass fraction of each species present in the field, . 
i.e. Hg, H, Og, 0, HgO, OH, Ng . 
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(V) 


turbulence energy and dissipation rate of turbulence 
energy. 

The program also predicts the distributions of auxiliary 
quantities; namely, temperature and density distributions, 
and drag and boat transfer at the walls. 

The program is written in standard FORTRAN IV language and, 
although it has been developed on a GDC machine, it 
will require only minor modifications to run on other 
machines with FORTRAN IV compilers . 

1 . 6 Outline of this Report 

The remainder of this report is divided into the following 
eleven Chapters : 

Chapter 2 is concerned with the mathematical formulation and 
physical models employed in SHIP. 

Chapter 3 then describes the finite-difference equations and 
outlines the numerical solution algorithm for solving the 
relevant equations . 

Chapters 4 to 8 are devoted to program topics and are intended 
primarily as the operational manual for the program. 

Emphasis is therefore placed on the description of the program 
structure, the function of the subrc^tines and on the input/ 
output features. These sections provide all necessary 
instructions to permit formulation of problems which fall 
within the scope of SHIP. It also provides information to 
aid interpretation of the predictions obtained from the 
program. 

In Chapter 9 ^ sample predictions for some test cases specified 
by NASA Langley are discussed. Suggestions for further 
extensions and refinements of SHIP are included in Chapter 10. 

Relevant literature references (Chapter 11 ) and nomenclature 
(Chapter 12 ) close this report. 

Appendix A describes in detail the equilibrium chemistry 
model used in conjunction with the computer program to predict 
the properties in a hydrogen-oxygen flame. 

Appendix B describes the Free Stream Boundary Conditions in 
supersonic flows . 

Finally, Appendix C gives a list of FORTRAN variables used in 
SHIP and a listing of the SHIP program is provided in Appendix 

D. 

The present report supersedes a previous CHAM report 
published as NASA CR-2655, March 1976 { 6 }. 
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2 


THE MATHEMATICAL AND PHYSICAL ANALYSIS 


2.1 Introduction 

This chapter details the mathematical and physical basis of 
the SHIP code. Section 2.2 is concerned with the system of 
coordinates chosen for use in the present work. Section 2.3 
is concerned with giving the differential equations for 
conservation of momentum, stagnation enthalpy and chemical 
species. 

Section 2.4 describes the turbulence model. Finally, Section 
2.5 provides the auxiliary information necessary to close the 
problem. 


2. 2 The Coordinate System 

The system of coordinates chosen for use in the present work, 
is quasi-orthogonal . The reason for this choice is that the 
representation of flows within domains whose cross-sections 
vary with axial position cannot conveniently be achieved 
through equations expressed in orthogonal coordinate systems. 
It is stipulated that whereas two of the coordinate axes 
(x,y) maintain mutual orthogonality throughout the flow- 
field, the third (z) is permitted to depart from orthogonality 
with respect to the other two, within specified limits. It is 
demonstrated below that these limits are consistent with the 
boundary-layer approximations {7}. 


The elements of the curvilinear system (5,n,C) are defined 
in terms of orthogonal, Cartesian coordinates (x,y,z) as 
follows: 


and 

Furthermore we define: 
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1 

1 

X X 

(a) 


y - Yg 

(b) 


^N ■ 


Z 

(c) 


= X 

w 

(a) 

ACj; 

= X “X 

e w 

(b) 

Arig 


(c) 

AHn 


(d) 


( 2 . 2 - 1 ) 


( 2 . 2 - 2 ) 
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The above definitions can best be appreciated with reference 
to Figure (4). The subscripts N, S, E and W, refer 
respectively to the North, South, East and West boundaries of 
the calculation domain in the x-y plane. 

The coordinates n and C are mutually orthogonal for all values 
of X,-. Furthermore, planes of constant 5 are approximated as 
planes of constant z. 


n 



FIGURE 4: ILLUSTEIATION OF THE QUASI-ORTHOGONiiL COORDINATE SYSTEM; THE 

CIRCLED CHARACIERS REPRESENT THE roRTH, SOUTH, EAST AND TOST 
BOUNDARIES OF THE CALCDLATICN DOMAIN 
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The above definitions for n,C make possible the use of 
coordinates that vary between 0 and 1 only. 

The components of velocity u, v and w are now defined as 
follows: u and v are normal to the y-z and z-x planes 

respectively, i.e. are aligned with the 5 and p coordinate 
directions, w is normal to constant C planes, but is 
permitted to depart from alignment with C by small angles; 
the limits of this inclination are prescribed below. The 
following mathematical consequences result from the above 
definitions. 


The coordinates (C,n,C) satisfy 


A » _L M + A 

3x 9C 9x 3n 


A = A M + A 

9y 9^ 9y 9ri 


A » A li + A 

9z 9^ 9z 9r) 


the general relationships; 


9n ^ 

9x 9c 9x 


In + A Ai 

9y 9C 9y 


(2.2-3) 


9 n 9 9 C 

9z 9? 9z 


which on application of the definitions of (2.2-1) and 
(2.2-2) reduce to: 


9 _ 1 A 

9x 95 


^ _1 9_ 

9y “ Anjj 9n 


(2.2-4) 


9z 





9 

95 


1 



9Ar)g ^ 9Arijg -.9 9 

9z 9z 9n 9C 


It can be deduced that, on applying the definitions (2.2-1) 
to relationships (2.2-3), the effect of non-orthogonality of 
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the C coordinate with respect to ^ and n is negligible only 
if the following conditions, are satisfied: 




(a) 


In . M 

dZ ? dz 


<< 1 


(b) 


(2.2-5) 


Conditions (b) are similar to the well-known boundary-layer 
assumptions; thus the definitions (2.2-1) along with 
conditions (2.2-5) permit the transport equations to be 
expressed in quasi-orthogonal coordinates, whilst retaining 
their boundary-layer character. 

The conditions (2.2-5(b)) are related to the area-ratio 
variation along the axis. It is obvious that this variation 
must be small to maintain the desired unstalled flow regime. 
Since stall, or axial flow recirculation violates the 
conditions required for boundary-layer flows (i.e. flows 
within the range of validity of SHIP), the above conditions 
are consistent with the physical nature of the flows under 
consideration . 

2. 3 The mean-flow Conservation Equations 

The differential equations listed in this section are the ones 
which express the conservation of momentum, mass, energy and 
chemical species in a three-dimensional flow with axial 
variations in cross-section. 

The form of these equations is restricted to flows which are 
classified as parabolic/hyperbolic. The term "parabolic 
flows" implies that: 

(a) there exists a predominant direction of flow (i.e. 
there is no reverse flow in that direction); 

(b) the diffusion of momentum, heat, mass etc. is negligible 
in that direction, and 

(c) the downstream pressure field has little influence on 
the upstream flow conditions. 

When these conditions are satisfied, the coordinate in the 
main- flow direction becomes a "one-way" coordinate: i.e. the 

upstream conditions can determine the downstream flow 
properties, but not vice versa. It is this convenient 
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behaviour of the parabolic flows that allows the use of a 
marching integration from an upstream station to a downstream 
one, for the solution of the following equations. 

a. Continuity 







(2.3-1) 

b. Transport of fluid property (b 



+ 

An^ an -^"dz 

+ 

£ 


Jr 9<J>1 + 1 L Jr Ml 

1 ♦ H J (An^)2 1 ♦ 5" J 

(2.3-2) 


In the above equation, (j> is taken to represent any fluid 
property transported by the flow, including the three components 
of momenta per unit mass, u, v and w. The qiiantities within 
square brackets are consequences of the curvilinear nature 
of the geometry confining the flow. The terms S, and T. 
represent, respectively the source (and/or sink) and the*^ 
diffusion coefficient for the transport of <t». 



The nature of the dependent variables and the associated forms 
of and F* are given in Table (1); and other symbols are 
defined in Section (12). The Table shows that the particular 
versions of (2.3-2) of interest here are: the momentum 
equations for each coordinate direction; the energy equation; 
an equation expressing conservation of total hydrogen; and 
the differential equations for turbulence properties. These 
equations together with the continuity equation (2.3-1) 
permit calculation of eight dependent variables, namely, 
u, V ,w,p , f , h ,k and e. 

The omissions from the equations are the shear stresses and 
diffusion fluxes acting in the z-direction. These omissions 
accord with the definition of parabolic flows and with the 
consequent necessity to ensure , that no influence from downstream 
can penetrate upstream. 

A further point to note, is that in subsonic flows the symbol p 
used for the pressure in the z-momentum equation is different 
from the symbol p in the two other momentum equations. This 
is a reminder of the fact that in the calculation procedure 
an inconsistency is deliberately introduced into the treatment 
of pressure, and that the quantities p and p are calculated 
differently. The pressure p can be thought of as a form of 
space-averaged pressure over a cross-section, and the gradient 
9p/3? is supposed to be known (or calculated) before 
calculation of the lateral gradients 9p/9n, 9p/9^. This 
practice is implicit in two-dimensional boundary-layer theories 
also. It is the final step to be made in preventing downstream 
influences from propagating upstream. If this step is omitted 
the solution is often wholly unrealistic physically. This 
inconsistency in the treatment of pressure, it may be said, 
is one part of the price to pay for making the equations 
parabolic; the gain is the freedom to employ marching 
integration and to use two-dimensional computer storage, even 
though the flow is three-dimensional and the full equations are 
elliptic. There is however ^ penalty for wholly supersonic 
flows, for which the three momentum equations can share the 
same pressure, without impairing the marching-integration 
feature . 

2 , 4 The Turbulence Model 

The effective turbulent-transport coefficients Peff* ^k > ^e> 

Fh and Ff are determined by means of a two-equation (k-e) 
model of turbulence {5}. According to this model, the 
magnitude of the viscosity depends only on the local values 
of the turbulence kinetic energy, k, on the dissipation rate 
of turbulence energy, e, and on the fluid density p. 

The form of the two transport equations for k and e is identical 
to (2.3-2); the transport coefficients and source terms are 
provided in Table (1) 
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Table (2), below, provides values for the laminar Prandtl, 
number, a, and for the turbulent Prandtl number 

Also given in Table (3) are values for the constants 
associated with the turbulence model: and C 2 which are 

required in the source term of the k and e equations; 
which is used to obtain the turbulent viscosity as will be 
described in Section (2.6) below, and K and E which are 
constants required in the law-of-the-wall formulation 
described in Section (3.4-4) below. 


TABLE (,2) 


Values of the laminar and turbulent Prandtl numbers 


<l> 


or 

k 

1.0 

.7 

e 

1.3 

. 7 

mm 

h 

. 9 

. 7 

f 

.9 

. 7 


TABLE (3) 

Values of the constants in the turbulence model 

formulation 


K 

E 

S) 

^1 

^2 

.42 

9.0 

.09 

1.44 

1.92 


Previous experience of applying the above equations to a 
large number of flow situations, has revealed that these 
constants tend to be nearly universal constants. However, 
a "fine-tuning" of their values for each particular flow 
situation may be necessary to optimise the results. 
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2.5 Auxiliary Information 

In addition to the partial differential equations, the 
complete specification of the mathematical problem requires 
provision of auxiliary information of four types: 

• Inlet conditions, i.e. initial values of dependent 
variables corresponding to the position of the coordinate 
along the predominant flow direction, at which solutions* 
to the set of equations are initiated; 

• Boundary conditions, i.e. conditions of all the dependent 
variables at the E, W, N and S boundaries, as a function 
of C »’ 

• Physical hypotheses which permit the calculation of 
diffusion coefficients as well as sources and sinks of 
each variable, in terms of the dependent variables of 
these equations, over the entire flow field; and 

• Certain relationships among the thermodynamic and 
transport properties. 

2i5-l Inlet conditions 

Information to start the marching integration is needed. 

This information must be provided as the inlet conditions of 
all dependent variables, fluid properties and the other 
auxiliary variables, at the plane at which solution is 
initiated. For the flows under consideration, this means 
specifying the velocities, pressures, enthalpies, species 
mass fractions, turbulence quantities and the physical 
properties, density, viscosity and specific heat. 

Any type of uniform or non-uniform distributions at inlet 
may be specified and supplied to the calculation procedure 
in a simple manner. For example, a distribution of 
experimentally-determined velocities may be supplied as a 
function of grid position. So also may the static-pressure 
and temperature distributions. 

2.5-2 Boundary conditions 

Boundary conditions along the N, S, E and W surfaces 

(see Figure 4) must be specified for each variable. Any one 

of the above boundaries may be either (1) a symmetry plane^ 

(2) a wall or (3) a free s\irf ace; and , as a consequence, 
different boundary conditions may apply on the different 
faces of the domain. These boundary conditions can be 
specified as the value of the variable <p or the flux of <j> 
through the surface. A detailed discussion of boundary 
conditions is given in Section 3.4.3. 
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2.5-3 Physical hypotheses 


• The gas density is calculated by the ideal-gas law, 
namely : 


; (2.5-1) 


weight W is calculated from: 


; (2.5-2) 


over all the chemical species 
j in the flow field. Details on the calculation and use 
of density in compressible-flow situations are given in 
Section 3.5. 

• The laminar viscosity is taken into account only in the 
vicinity of walls, through the "wall functions" {6} to be 
described in Section 3; 4-3. Its value is constant. 

There is no need for its inclusion in the core of such 
highly turbulent flows as the ones under study. 

• The turbulent viscosity is determined by means of the 

k-e model of turbulence employed. According to this model, 
the magnitude of the viscosity depends only on the local 
values of the turbulence kinetic energy k, on the 
dissipation rate of turbulence energy, e, and on the fluid 
density, p. The turbulent viscosity is then given by: 

= Cj^pk^/e ; (2.5-3) 

where Cj^ is a constant, given in Table 3. 

The length scale in this model is obtained from: 

k 3/2 

5- = Cjj I ^ , (2.5-4) 


RT 


where the mixture molecular 


1 » Z !!i 
* 3 


where the summation is taken 


and may be used for printing purposes, since its physical 
interpretation is more readily understandable than that 
of e. 
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2.5-4 Thermodynamic and transport relationships 

In order to calculate the fluid temperature and density, it 
is necessary to know the mass fraction of each chemical 
species at a given location. 

Thermodynamic equilibrium is assumed so that the species 
mass-fraction distribution is calculated from the local 
temperature, pressure, and element fractions. Details of 
this technique are given in Appendix A. 

Another important relationship is between the stagnation 
enthalpy and the temperature. This relationship is: 

2 . 2 2 

h = Z m.h° + ^ ■ + Cp (T - T^g^) + k (2.5-3) 


where Cp is the mixture specific heat and is the 

reference temperature for which the species enthalpy of 
formation, h. is defined. Strictly speaking the right hand 
side of equation (2.5-3) should include the local kinetic 
energy of turbulence as an additional term. This term has 
been neglected here because it is many orders of magnitude 
lower than the other three terms for the flows considered. 
For the present work, T „ is taken as 0°K so that C is 
defined by: ^ 


C H 
P 


T-T 


T h-h „ 

■ ^ “ T 

’’ref 


ref 


(2.5-4) 


where h is the specific enthalpy at temperature T. 
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3. THE NUMERICAL ANALYSIS 

3.1 Introduction 

In this chapter, details are provided of a nvunerical scheme 
reported by Patankar and Spalding {7} for the solution of the 
mathematical problem described in Chapter 2. The layout of 
the rest of this chapter is as follows: Section 2 describes 

the procedures adopted in the discretization of equations 
(2.3-1) and (2.3-2). The computational structure of the 
numerical scheme employed is outlined in Section 3. In 
Section 4, the manner of incorporation of auxiliary 
information into the computational procedure, is briefly 
outlined. The treatment of compressibility is reported in 
Section 5 which closes the chapter. 

3.2 The Discretization Procedure 

The finite-difference equivalents of the differential 
equations are obtained by integrating the latter over the 
control volumes which surround the nodes of a grid system. 

This procedure is described in the following sections. 

3.2- 1 The Grid system 

The finite-difference grid used consists of; 

(a) In the n-C planes, a system of intersecting, 
orthogonal grid lines of constant n and C . No 
restrictions are placed upon the spacing between the 
lines in any given direction. 

(b) Planes of constant C at which solutions are obtained 
are arrived at by taking succesive increments (i.e. 
forward steps) along the C direction (main flow 
direction). 

The limits on the size of the forward step, are 
governed by considerations of stability and accuracy 
of the numerical procedure. 

3.2- 2 Storage locations 

The intersections of the grid lines mentioned above are 
termed grid nodes. All the fluid properties with the 
exception of the velocity components u and v, are stored at 
the grid nodes. The velocity v is located midway between 
grid nodes in the n-direction, and velpcity u similarly 
located along the C direction. Figure (5) illustrates this 
"staggered" grid system in the n-C plane. The boomerang- 
shaped envelopes shown on the figure enclose the triads of 
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FIGURE -5; STAGGERED LOCATION OF VARIABLES 



points denoted by the single letter, N,S,E,W or P, 
which represent a unique computer storage location.* 

In the algebra connected with the discretization procedure, 
when a variable is required at points other than those at 
which it is located, averages of neighbouring values are 
used to arrive at the value of the variable at that point. 

3.2-3 Control volumes 

The control volume surrounding each grid node P, indicated 
• in Figure {6} has two faces that coincide with two cohstant- 
t, planes. One of these, at which integrations of the 
partial differential equations are to be performed, is 
designated the downstream (D) station. The other faces are 
midway between the nodes, so that the velocity components 
giving rise to convective fluxes along the n and ^ 
directions are located bn the faces themselves. 

Figure {?}, illustrates this point, with reference to a 
cross-section of the control volume in the g-h plane. It 
is over such control volumes, that balances of w, (J> and 
mass are made in the calculation procedure. Similar 
control-volumes , resulting from the "staggering" of locations 
on the grid, are defined to surround the locations of the 
V and u velocity components. Three sets of control volumes 
are thus identifiable over the entire calculation domain. 

A slight modification to the variable location and 
control -volume definitions is made in the region of 
boundaries of the calculation domain. The cpntrol 
volumes corresponding to the near-boundary velocities v in 
the case of N and S ,and u in the case of E and W boundaries , 
are arranged to extend right up to the boundary. Figure (8) 
illustrates this point. 


* Note on change in notation. N,S,E,W were used in 

Chapter 2 to denote domain boundaries. From here on they 
will denote neighbouring grid nodes. 
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FIGURE 7 ; CONTROL VOLUME WHEN THE GRID EXPANDS TO 
ACCOMMODATE EXACTLY THE DIFFUSER SHAPE. 
NOTE THE INCLINATION OF THE w-VELOCITY 
ARROW WITH GRID LINES 
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near - boundary 
V- control volume 


FIGURE 8; THE NEAR-BOUNDARY MODIFICATION TO CONTROL VOLUMES 
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3.2-4 The Discretized equation 


The general discretized equation is obtained by 
Integrating equation (2.3-2), for each variable, over 
the appropriate control volume. In the integration, 
the following assumptions are made about the manner in 
which the variables vary between grid nodes. 

(a) In the z-direction, <j> varies in a stepwise manner; 

i.e. the downstream values of (j> are supposed to prevail 
over the interval from the upstream station to the 
downstream one except at the upstream station itself. 
This makes the present finite-difference scheme a 
fully-implicit one. 

(b) For the calculation of the z-direction convection, and 
of source terms that may depend on <j) , the variation 
of <j) in the X 5 ^ plane is also taken to be stepwise. 

Thus, in the xy-plane the value of <}> is assumed to 
remain uniform and equal to <{ip over the ■ control volume 
surrounding the point P and to change sharply to 
(j)jj,<J)g, (J)g , or (fiyj outside this control volume. 

( c) For the cross-stream convection from the yz and xz 
faces of the control volume, the value of (f) convected 
is taken to be the arithmetic mean of the (j) values 

on either side of that face, except when this practice 
is altered by the "high-lateral-flux" modification 
mentioned below. 

(d) For diffusion across the yz and xz faces of the control 
volume, we assume that cp varies linearly between grid 
points, except when the high-lateral-flux modification 
dictates otherwise. 

The result of these operations is an algebraic equation 
for each grid location, representing the discretized 
form of the balance of the variable, over the control 
volume corresponding to that location. For a general 
dependent variable <j>, .this equation takes the form: 

*P.U * <*N " *P> - *s^2<»P " 

* fe (1‘E * <*P * ■* Sp* 

* ■'S - '♦p^ - ’'s <*P - ♦s’ 

* (*E - ♦p) - <*P - fvj> (3.2-1) 
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The factors f^, f , and f^ are interpolation factors 
to account for the staggered cells. For the <j> cells, 
they are all 0.5. 

Attention is drawn here to the manner of discretization of 
the convective terms in the n- and f directions. The 
method used is designated the "high-lateral-flux 
modification" and is discussed in detail in {8}. It 
is briefly described later. 

The coefficients of equations (3.2-1) are now defined as 
follows : 


a 



(a) 


" [pv -pw. (Aiig) * 

A5pACA5g 

(b) 

L? 

1 


AripA ^ AtIj^ 

Ji 

(c) 

LP 

= l4 n “ - l’^ + 

nP,U e w n s 


(d) 

t5 

AnpAcAn*! 
* r. 


(e) 





<t>,J 6rijAr)jj 


(f) 


= S^^pA?pAnpACA?gAri^ (g) 

where , 

A's and 6's represent the widths of control-volume faces 

and internodal distances respectively (see Figure 9)» 

and i and j stand for locations (e,w) and (n,s) respectively. 
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Ac is the forward step size, and is the integrated 
form of the source (and/or sink) of <j> , always expressed 
for reasons of numerical stability in a linearised form 
as : 

Sp^ = Sy + Sp <J)p (3.2-3) 

Where Sy and Sp are calculated from values already in 
computer store at that time. Expressing source terms in 
linearised form enhances the stability of the finite- 
difference equations, and is recommended as a useful 
practice to be applied wherever possible. Further 
it is alsc necessary that Sp be negative always if 
stability is to be ensured. 

On re-arranging the terms of equation (3.2-1), one 
obtains; 


S A.’ 
1 


<^1 ^ 


(3.2-4) 


i = N,S,E,W,U i = N,S,E,W 

where the coefficients are defined as follows :- 


II 

< 



(a) 

II 

< 

* 1 


(b) 

A ' = 

m tl 1 T ^ 

^n 2 '"n 


(c) 

A ' s: 

^ 1 "s^ 


(d) 

A ' = 



(e) 

CO 

II 

(Sy + Sp(j)p) 


(f) 


where the L’s represent convective contributions, T's 
the diffusive and S the source (and/or sink) contributions 
to the balance of (j>. 

The subscripts are associated with points on the grid 
system and superscripts to the co-ordinate directions 
for which the coefficient is appropriate. The location 
of points e,w,n,s, E,W,N,S is given in Figure (9). 
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Note that in the derivation of , account has been taken of 
the fact that the faces of the main control volumes pass 
midway between grid nodes so that the interpolation factors 
for the calculation of variables on these faces are equal to 
0.5. This is not true for the control volumes appropriate to 
u- and V- velocity components in the general case of non- 
uniform spacing between grid nodes. Direct interpolation 
is then necessary so that the value of u at point P is given 
by : 


“e * “p 


and similarly for V. Note also that: 


w 


= 1 - f 


i 

It is possible for the convective contribution Lj of the 
coefficient A' to become large on occasion, resulting in 
the coefficient becoming negative and causing physically 
implausible results when equation (3.2-4) is solved with 
such coefficients. The high-lateral-flux-modification, 
mentioned above, is introduced to overcome this possibility. 
This scheme consists of replacing all the coefficients of 
the form T^by T^ as follows: 




_1 

2 






(3.2-6) 


where | A | signifies the modulus of A. 


3.2-5 Equations for velocity components 

The finite-difference equations for the velocity components 
are described by the same form as (3.2-4) but contain an 
additional source term representing the pressure-gradient. 
Since we shall be using the velocity equations later in 
deriving the continuity balance, it is necessary to note 
here their form. 


5-direction momentum equation 

t-aW . ciW , 9P 

Wp = EA. w. + Sp + Dp 3^ 

i=E,W,N,S 

n-direction momentum equation 

P = 2AYv. + ^ d; (Pp 

i=E,W,N,S 


- Ps) 


(3.2-7) 


(3.2-8) 
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(3.2-9) 


g-direction momentum equation 

Up = + Sp + Dp ” ^W) 

i=E,W,N,S 

where, 

denotes again coefficients of the form: 


At 

1 

ea!^ 


(3.2-10) 


i=E,W,N,S,U 



where 0 

now stands for u, v 

and w* 



= -A^AnpA^pArij^ACg 

, (a) 



= -AcACpACj; 

, (b) 

(3.2t11) 


= -ACAnpApj^ 

, (c) 



and, 

Sp represents the source (and/or sink) of each 

velocity excluding the pressure-gradient terms. 

The pressure term employed in the w-equation depends on 
whether the local flow is subsonic or supersonic. When 
the flow is subsonic, it is necessary , '{7} , in order to 
render the equations parabolic, to write the source term 
as a mean pressure gradient, i.e., 

= If = ^ ~ % (3.2-12) 


Note that these A.*s now relate to the staggered control 
volumes described^earlier . 



where p denotes the mean pressure which is determined for 
a confined flow from the_requirements of overall continuity; 
for an unconfined flow, p is simnly the free-stream pressure. 
For a supersonic flow in contrast, it is possible, without 
destroying the parabolic nature of the equations, to employ 
the local pressure to calculate the pressure gradient, i.e., 

S = -^ = ^2jIL (3.2-13) 

dz Sz 


This practice, it should be observed, allows full account 
to be taken of pressure waves in supersonic flow. 


3.2-6 The continuity equation 

The derivation of the finite-difference form of the 
continuity equation is quite simple. It merely states 
the requirement that the inflows and outflows of mass 
are locally in balance at all the grid nodes in the flow 
domain. In the present solution procedure, the requirement 
of mass continuity is satisfied by correcting the pressure 
field via a pressure-correction equation. The details of the 
latter will be given in Section 3.3 below; here we state 
merely the requirement of continuity as follows: 

continuity 


iol - G^) AnpAt An^ {gJ - G^) A?pA5 A5^- 

(3.2-14) 


where , 


represents the mass velocity along direction j at location 
i. For example: 



pw' { 


dz 



Ahp y and A^p y represent the 
ACp respectively. 




Js 

upstream values of An 


(3.2-15) 
p and 
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It should be noted here that two distinct densities are 
computed and used, at each step. The downstream density 
(Pp) is used only to compute the mass flux in the ? 
direction at the downstream face of the cell (thus it 
appears only in the last term of equation (3.2-14); 
the upstream cfensity ( Pp y) is used for calculating lateral 
mass fluxes, and the upstream axial mass flux. The two 
densities are calculated as follows: 


P 


P 


Pi,W. 


RT 


P P.U 


P,U 


(3.2-16) 




Pp ,U ,UU ' 


RT 


P, UU 


(3.2-17) 


where Wp is the local mixture molecular weight, R is the 
universal gas constant, and subscript 'UU' refers to the 
plane two steps upstream of the plane in question. 

3.3 The Solution Procedure 

3.3-1 The computational algorithm 

The above set of finite-difference equations has to be solved 
for all the variables simultaneously, at the downstream 
station. After completion of the solution at the downstream 
station, a forward step in the ^-direction is taken 
and the procedure repeated. 

The numerical algorithm embodied in SHIP is the SIMPLE 
(for Semi-_Irap licit Method for Pressure-M nked Equations) 
scheme reported earlier in {7} by Patankar and Spalding. 
Central to this scheme is the idea of seeking a non-iterative 
marching-integration procedure that takes full advantage of 
the boundary -layer character of the flow field. 

At each axial location the variables are computed solely 
from values at the upstream location; no reference is made 
to downstream properties. Further, the coefficients (A's) 
used in the equations are evaluated from the values currently 
in store; for the velocities , this means use of the upstream 
values. The sequence of calculation steps to evaluate the 
flow properties at any axial station is as follows: 
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1) First the pressure field p and the mean pressure p 
at the axial station considered are assigned 'guessed' 
values. The general practice is to employ the 
calculated upstream pressures as the guessed values 
an,d to estimate p using the upstream pressure gradient 
dp/dz. 

2) The three momentum-equations are solved to get a first 
approximation to the velocity field at the longitudinal 
station. 


3) 


The resulting "starred" velocity field is used in 
conjunction with the discretized continuity equation 
to arrive at a distribution of "pressure-corrections" 
p' as follows.* 

(a) First the pressure and velocity fields are 
expressed as : 

P = p* + P' 


u = u* + u' 


(3.3-1) 


V = V* + v' 

where the primed quantities represent the 
corrections to the approximate *starred* values. 
The latter will not, in general, satisfy the 
continuity equation, but will give rise to a net 
mass source at P. 


(b) 


It is now required to obtain the corrections to 
the velocities and pressures so as to reduce this 
mass source to zero. To this end, the substitution 
of equation (3.3-rl) into the momentum equations 
results in the following expressions: 


"P “ ^P ^Pp " Pw> 

= Sa];v’ . d; (p- - P-) 


(3.3-2) 


where the summation Z is carried out over the grid nodes 
neighbouring the corresponding velocity. 


♦Readers not interested immediately in the details of the 
pressure-correction equation may skip to step 4, but should 
note that the pressure- correction equation has the same 
general structure as that of (3.2-4), and that it is solved 
in the same manner as the others. 
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(C) 


The relations expressed above are substituted 
in the finite-difference form of the continuity, 
equation (3.2-14); and the coefficients for pp 
are collected and rearranged. The pressure 
correction equation has then the following form: 


Pp 





(3.3-3) 


i=E,W,N,S 


lit I 

where Pg, p™, Pj^ and Pg are corrections to the 
pressures at the nodes E,W,N and S, the A.'s 
involve p*s, D's and other geometric quantities 
appearing in equation (3.2-14), and the mass 
source.i (or mass Imbalance) has been incorporated 
into Sp 

For compressible flow, account must be taken of thi 
effect of a pressure change at P on the mass flux 
at the downstream face of the cell. Thus for 
supersonic flow the mass-flux change ( pw) p 
related to the pressure correction Pp by: 

(pw)^ = {p*pD^ + Pp (3.3-4) 


where, Pp and Wp are computed from the guessed 
pressure Pp; (dp/dp)p is, from equation (3.2-16) 



W. 




RT 


P,U 


(3.3-5) 


and Dp 



(3.3-6) 


and is deduced from the finite-difference form of 
the momentum equations. 

For subsonic flow, equation (3.3-4) must be modified 
to account for the fact that Wp no longer depends 
on Pp; the modified expression is: 


(pw) 


P 


(^) 

'•dp^P 


'P 


Wt> Pp 


(3.3-7) 



(d) Equation (3.3-3), with the coefficients modified 
to account for the above-descrihed influences of 
Pp on (pw)p, is solved in the same manner as the 
momentum equations. The resulting pressure- 
corrections are then used in correcting the 
pressure and velocity fields simultaneously. 

4) The equations for the remaining variables (i.e. enthalpy, 
turbulent energy and its dissipation rate, and mass 
fraction of hydrogen in any form) are then solved, by 
using the corrected velocity fields. 

Steps 1-4 complete the operations at a given downstream 
station. A new step is then taken and the process 
repeated, until the region of interest is covered. 

3.3-2 The solution of the equations 

The linear algebraic equations are solved by performing 
repetitive sweeus of the tri-diagonal matrix algorithm 
(TDMA) along the two cross-stream directions n and 

The equations are solved along lines of constant p and 
of constant and in doing so, the variables located at 
adjacent lines are kept constant. Thus for the p-direction 
sweep we have: 

<Pp = (3.3-8) 

where the terms enclosed in the parentheses are assumed to 
be fixed. 3imilarly, for the ^-sweep we have: 

<j>p — ^3*^3 (3.3—9) 

where now the values and (bg are kept fixed. The 

terms contained in 3^ however, are always like the A coef- 
ficients, evaluated Trom variables in stoi'e before the TDMA sweep 

The number of TDMA sweeps required to obtain an accurate 
solution to the equations depends on the equation being 
solved. It has been found from ' experience that for the 
pressure-correction equation it is necessary to perform 
more sweeps than for the rest. The reason for this 
increase is as follows. The equation for a general variable 
(t> is dominated strongly by the contribution from upstream 
(i.e. the term u 4' p y in equation (3.2-5 (f)). Therefore 
the coef ficients^Ajvf etc, are relatively small, and use 
of somewhat approxxm^e values for 4*5 etc, does not 

result in significant errors. For the pressure-correction 
equation however, there is no contribution from upstream, 
and the pressure correction at P is related strongly to 
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the N,S,E and W values. It is therefore necessary to 
per f oral several sweeps to obtain accurate values,' 

3.4 Boundary Conditions 
3,4-1 The solution domain 


SHIP is set up to make calculations in a rectangular-sectioned 
domain of any aspect ratio. The boundaries of this domain may* 
be either solid walls, symmetry planes or free boundaries; 
and indeed any combination of the above boundaries can be 
handled. Furthermore, for walls, the distance of each wall 
from a reference plane can be specified as an arbitrary, 
but smooth, function of distance along the principal flow 
direction. 

In the region of f ree-boundaries the flow is unconfined and 
expands as it oroceeds downstream. Two considerations must, 
therefore, be taken into account here, namely the prescription 
of the rate of expansion of the solution domain and the 
specification of the mass velocities at the outer free- 
boundary, The latter will be described in Section 3.4-3 below. 

For the former, a slope of the boundary is prescribed; and if 
this leads to solution showing excessive velocity and 
temperature gradients the calculation is repeated by moving 
the f ree-boundaries outwards or changing the slope of the 
free surface in the desired direction to diminish gradients. 

3.4-2 General policy of treating boundaries 

The SHIP program requires the specification of boundary 
conditions on the N, S, E and W boundaries to ? - q planes 
of the calculation domain. 

A clear distinction is made, for all dependent variables, 
between boundary values and values internal to the domain. 

The main machinery of the program leaves the boundary 
values unchanged, although it uses them in determining the 
internal values. Thus, the procedure is so structured that 
it nominally solves the fixed-boundary-value problem. When 
boundary values are not known however, appropriate modi-j- 
fi cations are devised which permit the single structure to 
be used. The following sections describe such modifications. 

In general, boundary-condition information can be supplied 
to the numerical calculation procedure in one of four ways. 

The boundary values of the dependent variables themselves 
can be modified; or the values of T at the boundary nodes 
can be suitably adjusted. Alternatively, the source terms 
for the near-boundary control volumes or the finite-difference 
coefficients themselves can be suitably modified, SHIP is 
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equipped with source-term modification practices. 

3.4-3 Treatment of the types of boundaries handled by SHIP 

SHIP allows for any of the four boundaries to be (1) a 

symmetry plane (2) a wall, or (3) a free surface. 

(1) Treatment of symmetry planes 

At symmetry planes and at the axis of symmetry, the 
velocities normal to the boundary are zero; and there 
are no fluxes across them of other flow variables. In 
the present solution scheme such boundary conditions 
are incorporated as follows. The normal velocities 
at these boundaries are prescribed to be zero a priori 
and not altered; for other variables, because both the 
convective and diffusive fluxes are zero, the 
appropriate finite-difference coefficients connecting 
the boundary node to the near-boundary one are set to 
zero, thus breaking their links. 

(2) Treatment of wall boundaries 

a) At solid walls, all the three components of velocity 
are zero and are prescribed a priori. For turbulent 
kinetic-energy, k, and its dissipation rate, g, the 
boundary-conditions are prescribed through wall 
functions described below; no reference is made to the 
values of k and e on the wall nodes. The values of k 
and e on the wall therefore, are prescribed arbitrarily 
to be zero, and have no^ influence on the solution scheme. 
For enthalpy, the SHIP code provides the choice between 
specification of an adiabatic wall or a constant 
temperature one, for each wall. In the former case, 
(prescribed zero-flux boundary condition)^ the value 
stored as wall temperature does not enter the 
calculations . 

In the region close to the wall, the correct fluxes of 
momentum are calculated through wall functions described 
below. 

b) Wall functions 

The expressions for T appropriate to turbulent flow 
are not strictly valid in the vicinity of wall 
boundaries to the flow, where laminar viscosity plays 
an important role. If the near-boundary grid nodes 
are sufficiently far away from walls, the turbulent 
viscosity can continue to be used for internal grid 
nodes. However, means must be provided for 
calculation of the correct shear stresses as well as 
fluxes of other dependent variables at the wall 
boundaries. 
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Provisions for such calculations are made in the program, 
and use the so-called wall-function concept. In this 
concept, the flux of a variable «() at a wall boundary, 
is expressed as : 

AM = r ^near wall - ’^wall (3.4-1) 

'’“wall ^(f>,wall '^wall 

where denotes the normal distance from the wall 

to the near-wall point. Values of T. ,, are obtained 
from the presumption that in the regldn adjacent to 
wall boundaries, the dependent variables obey, for 
turbulent flow, a modified form of the semi-logarithmic 
law-of-the-wall. The. formulae used to calculate 
F(^ wall provided in Table (4). 


TABLE 4. Diffusion Coefficient Formulae 



^(j), wall 

Velocity components 
normal to the wall 

0 

Velocity components 
parallel to the wall 

> 11.5 : 

^ J-n {Ey*^} 

< 11.5 ; U 

K 

0 

e 

— 

All other <j>'s 

( only for h 
in SHIP) 

y".11.5:.^ '-i 

°t,«. «n {Ey ) + P ) 

< 11.5: ^ 
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The definition of y"*^, in the table, is a 
generalization due to Spalding {8} of the conventional 
form, in that; 


+ 

y 


,pCp^k^6 

ii 


(3.4-2) 


where 6 denotes the distance from the wall, at a 
location with which the values of p and k are 
associated. The constants K and E are obtained from 
the conventional form of the law-of-the-wall : 

u'*’ = ^tn {Ey'*’} (3.4-3) 


and are given values of 0.42 and 9.0 respectively. 

The boundary conditions of k and e are provided as 
follows : 


The diffusion of kinetic energy k, to the wall is known 
to be negligible and is set to zero and a balance 
equation for k, regular in other respects, is Solved for 
control volumes adjacent to wall boundaries. The 
diffusion of dissipation rate e to such boundaries is 
more difficult to express. Instead of the attempt to 
calculate F , , , use is made of the knowledge that 
the length Sd^re'*'i2< varies linearly with distance from 
the wall, in the neighbourhood of the wall. The 
dissipation rate is then calculated from this length 
scale from; 


e = C ^ ^ (3 4-4) 

^near wall 4 ) K6 

The practice adopted is to fix e wall above 

value, without disturbing the general"caIculation 
procedure, in the following manner; 

^ ®near wall (3.4-5) 

-L (b) 

30 

is a large number, say 10 



where L 
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The expression for F. ™qii ^ general dependent 

variable (f) is based* An xne expression indicated in 
Table (4), the value of being calculated from: 


P. = 9.24 

<p 


where Oj. and a. 


t,(j) 


- 1 



• * 



the 

laminar 




(3.4-6) 


values ^f Pran^tlf /Sdhmidt number appropriate to the 
transport of (}) . The above equation is used in SHIP 
only for heat transfer. 


3. Treatment of free-stream boundaries in supersonic flow 


When the free-stream is supersonic, the velocities 
normal to a free-stream boundary can be deduced from the 
following equation: 


V 


D-D. 


p w 
*^00 00 



(3.4-7) 


where: w^^/YPco (3.4-8) 

and «> denotes undisturbed free-stream conditions. 

The above equation is based on the assumption that all 
changes to the lateral velocities, from their free-stream 
value along the free-boundary, have been brought about by 
a succession of small waves emanating from the viscous 
region. It has been programmed to compute the lateral 
boundary velocities, normal to the free boundary under 
consideration (i-.e. the v-component for the North and 
South boundaries, and the u-component for the East and 
West ones-. 

Note that: 

1. The formula is valid only for > 1 

2. The formula assumes v^ and u^ equal to zero 

3. In the computer program, p in (3.4-7) is actually, 
p . , where nb denotes values at the £ear-boundary 
grid point. 



ir 


4, The density at the boundary node is calculated 
from but otherwise from conditions at «>. 

The formula holds true whether fluid is caused by the 
pressure changes to leave or enter the free boundary. 

If fluid enters, its stagnation enthalpy should be that 
of the free stream. 

The longitudinal boundary velocity is calculated from the 
following equation derived from isentropic considerations: 


(3.4-9) 


where the subscript b denotes values at the boundary. 

The derivation of equations (3.4-7) and (3.4-9) is given 
in Appendix B. 
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A DESCRIPTION OF THE COMPUTER PROGRAM 


4.1 Introduction 

In this chapter an overview of SHIP is given with a 
description of the program flow and control, and 
definitions of important variables. 

4.2 General Structure 


The general structure of the SHIP computer program is 
shown in .Figure (10). The divided, vertical box on the 
left is the main program, MAIN, which has the chief function 
of controlling the calling sequence of the subprograms. 

The starting point for SHIP is at the top of this box. 

Points where calls of subprograms are made in SHIP are 
denoted by horizontal lines. Arrows indicate the "flow" 
of the computer program to and from the subprograms which 
are shown in boxes. The major calculations and loe^ical 
decisions made by SHIP are indicated in the spaces between 
the horizontal lines marking calling points in SHIP. 

The subprograms are contained in subroutines which are not 
individually called; and all start with an ENTRY statement. 

The subroutines are therefore never "called" by their 
names, but by the name of the appropriate member subprogram. 
(Thus CALL AUX is meaningless; but CALL SOURCE is a correct 
statement). No subprogram has any argument; the information 
is everywhere transferred through COMMON. 

The SHIP program consists of 8 subroutines; they are given 
the names BLOCK DATA, MAIN, ALLMOD, AUX, PRINT, SOLVE, STRIDA 
and STRIDE. These subroutines can be classified in three 
categories: problem-dependent, physical-modelling-dependent 
and invariant. BLOCK DATA and ALLMOD are the problem-dependent 
routines, i.e. they provide for specification of inlet 
conditions, boundary conditions, geometry, fluid properties, 
etc. AUX forms the physical-modelling section of the program; 
in this the auxiliary quantities such as density and 
viscosity are calculated from physical laws and from the 
turbulence model. SOLVE, STRIDA, STRIDE and PRINT are the 
invariant portions of the program; in SOLVE, STRIDA and 
STRIDE, the calculation steps of section 3 are programmed, 
and in PRINT, instructions are provided to print out the 
distributions of the flow variables. 
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START SHIP 


I NO 


CALC- 

ULATE 

MASS 

SOURCES 

AMO 

STEP 

LENGTH 


CALC- 

ULATE 

OVERALL 

BALANCES 


ISTEP 

NPJUMP? 

INSTEP = 
IINJ OR 

MODIFY 
DZ AND 
NSWP 
ISTEP= 
LASTEP?] 


IS 

GEOMETRY 

CHANGING? 

ISrEP^O? 
ITERATE 
*e) U, V, w| 


ZU> 

ZLAST? 



• Call STRID3 once 

• Call STRID4 for each 
dependent variable 


SAME AS STRID4 ABOVE 


FIGURE 10: GENERAL STRUCTURE OF SHIP COMPUTER PROGRAM 
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The subprogram MAIN does not fit neatly into any of these 
categories. It has a general (i.e. problem dependent) 
structure, but care is necessary to ensure that its 
general structure is kept intact. The routine MAIN first 
calls routines STRIDO, STRIDl and BEGIN. Through these 
calls', the initial information of geometry, grid and inlet 
profiles is prescribed. A loop for forward steps is then 
started. In this, calls are made in sequence, to routines 
UPSTRM, DENSTY, VISCOS, GEOMOD, STRID3 and STRID4. By 
these calls, the variables u,v,w,k,e and h at one down- 
stream location are calculated. STRID3 and STRID4, in turn, 
call routines dealing with the physical modelling section 
AUX and the problem-oriented subroutines in ALLMOD. Yfhen 
the finite-difference equations are assembled, through calls 
to AUX and ALLMOD, SOL'VE is called to perform the TDMA sweeps. 
STRID3 contains four sections each having identical calling 
sequences of routines placed under AUX and ALLMOD. STRID4 
calls the same routines as called by individual sections of 
STRID3; and solves the equations for the scalar variables 
k, e and h. 


The remaining part of this chapter gives a general overview 
of the computer program while succeeding chapters detail the 
operation of each subroutine. 

4.3 Some Programming Conventions and FORTRAN Equivalents of 
Main Variables 


A few of the FORTRAN names used in SHIP, which will be required 
in the following sections, are introduced here. The deoendent 
variables <t> are stored in an array F, which it is convenient 
to consider as a three-dimensional array F (I,J,NV). Here I 
and J denote the location (respectively along the 5 and n 
directions) and NV identifies a particular variable. The 
three velocity components and the pressure-correction are 
included in the F array; however, for ease in understanding, 
separate arrays U,V,W and PP are also used and made equivalent 
to parts of F as follows: 


U(I,J) 
V(I,J) 
W(I,J) 
PP(I , J) 


F(I,J,NVU) 
F(I,J,NVV) 
F(I, J,NVW) 
F(I,J,NPP) 


(4.3-1) 


The identifiers NVU, NVV, NVW and NPP, also used to identify 
U,V,W and PP elsewhere in the program, are assigned values 
2,3,4 and 1. The largest values of I,J and NV for which 
storage is provided in the program are denoted by IMAX, JMAX 
and NNV respectively and assigned values 12, 20 and 9. 
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Although two- and three- dimensional arrays have been 
mentioned above, the computer program formally uses one- 
dimensional arrays, v/hose subscripts are calculated each 
time they are required. Thus, F(I,J,NV) is referred to 
as F(IJNV), where: 


IJNV = I + JM(J) + NFM(NV) (4.3-2) 

the arrays JM and NFM being calculated once and for all 
from: 

JM(J) = (J - 1) * IMAX 

(4.3-3) 

NFM(NV) = (NV - 1) * IMAX+JMAX 

Also W(I,J) is referred to as W(IJ) where 

IJ = I + JM(J) (4.3-4) 

The four neighbouring points of the location IJ, corresponding 
to the points of the compass are referred to as UN, IJ^, 

IJE and UW. These are calculated as: 

UN = IJ + 

US = IJ - 
IJE = IJ + 

IJW = IJ - 

It is easy to see that points further removed become: 

IJNE = UN + 1 
USW * US - 1 etc. 


IMAX 

IMAX 

1 

1 


(4.3-5) 


For a given problem, all members of the F array for which 
provision is available, may not be required to be solved 
for. Furthermore' some of these variables may be obtained 
from algebraic equations and not from solutions of the 
partial differential equations. To provide for these 
alternatives, use is made of an array ISOLVE(NV). For 
ISOLVE(NV) equal to zero, the differential equation is 
not solved; solution is obtained for values of ISOLVE(NV) 
greater than zero. It is left to the user to make further 
use of this facility. 


Other arrays directly related to members of the F array are: 
IPRINT(NV), TITLE(K,NV), FLUXN(I,NV), FLUXS(I,NV), FLUXE(J,NV), 
FLUXW(J,NV). Values of F(UNV) are printed out if IPRINT(NV) 
is equal to 1; otherwise a printout is not obtained. 
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FLUXN(I,NV), stores the fluxes of the variable NV, on 
the North boundary; the others correspond respectively 
to the South, East and West boundaries. TITLE(K,NV) with 
the integer K taking values 1-9, stores a 36 character 
alphanumeric title of the corresponding variable NV. This 
is used to identify the values of the variables in the 
printout. 

The quantities k, £, h, f, and T form the remaining 
variables of the F family. They are identified by the 
indices NVK, NVD, NVH, NVF and NVT respectively and occupy 
NV locations 5, 6, 7, 8 and 9 in the array F(IJ,NV). No 
equation is solved for T; however, it is convenient to 
store it at F(IJ,9). Quantities p, p and T are stored as 
P(IJ), RHO(IJ) and GAMCIJ). They are not members of the 
F array, but are included along with F in a COMMON 
statement, and immediately following it. They can therefore 
be referred to as the (NNV + l)th^ (NNV + 2 )th and (NNV + 3 )th 
members of the F array, as Indeed they are in subroutine 
PRINT. The integers NVP, NRO and NMU are used to identify 
them. 

Although, in general, the values of T are different for 
each dependent variable, storage provision for only one 
set of r's is retained. Hence, at a given stage in the 
program, the GAM array contains the values of F, appropriate 
to the dependent variable under consideration. 

Similarly, provision is made, for the retention of one 
set of coefficients of the finite-difference equation at a 
given stage.. Thus, AXP(IJ), AXM(IJ), AYP(IJ) AYM(IJ) and 
Az7Ij) respectively represent the coefficients aI, A', A', 

Ag andAy in equation (3.2-4). " 

Similarly, SU(IJ) and SP(IJ) represent Sy and Sp respectively 
in equation (3.2-3). The quantities like dH in'^equation 
(3.2-9) are stored as DU(IJ), there being similar arrays 
DV(IJ) and DW(IJ) for the corresponding quantities 
associated with the v- and w- momentum equations, 

4.4 Program Control 

Program Control, including start, internal monitoring and 
stOD functions, is achieved through subroutine MAIN. 

Using the information supplied through BLOCK DATA, the first 
part of MAIN calculates and assembles all the geometrical 
information about the grid system through calls to STRIDO 
and STRIDl. Following this, in the same part, a call to 
BEGIN supplies the initial conditions for all the 
dependent variables. Having provided the starting conditions. 
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MAIN then enters the marching-integration loop. 

This part begins with the statement: 60 CONTINUE, 

Next, mass sources CSUM), step length (DZ) 
and the location of the station for the new calculation 
(ZD), are determined. This is followed by calls to SPECIE, 
DENSTY and UPSTRM. SPECIE calculates the chemical specie 
composition from element compositions, temperature and 
pressure. DENSTY is another member of the physical-modelling 
subroutine AUX, wherein the fluid density p, is calculated. 
UPSTRM stores the upstream values of the axial mass-flow- 
rate (GZ(IJ)), of the pressure (PU(IJ)) and of the height 
and width of the calculation domain (BYNU and BXEU) required 
by the calculation procedure. 

After calculating and printing the fluxes of the major F 
arrays to check their overall balances, a call is made to 
VISCOS which calculates the diffusion coefficients for 
momentum and modifies them properly to account for the 
presence of solid walls if necessary. 

If ISTEP equals NPJUMP, profiles of variables specified 
by IPRINT(NV)=1 will be printed out through a call to print. 
Otherwise, only overall balance information and station 
location will be printed out. 

The next section is for the calculation of the velocity and 
pressure fields. First a check is made to see whether 
injection through either or both North and South boundaries 
occurs. If this test is positive, a call is made to INJMOD 
(for the South boundary) and/or INJMOT (for the North 
boundary) which provide the appropriate boundary conditions 
through source terms. It should be mentioned here that 
no restriction is placed upon the relative simultaneous 
spacing of the North and South injectors. Next a check is 
made to determine whether the station for calculation is 
immediately downstream of injection. If this test is 
positive the number of sweeps on the pressure correction 
equation and' the step length are modified to increase 
accuracy and stability of the program. The calculation is 
stopped at this point if ISTEP > LASTEP. 

Next a check is made to determine whether the flov; domain 
under consideration presents varying geometry with axial 
distance or not. If this is true (IARCH=2), a call is made 
to subroutine GEOMOD and STRIDl, in order for the geometrical 
changes to be taken properly into account. If the domain 
of interest retains constant cross-section along the axial 
distance, the above two subroutines need not be called. 
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After this test, another is performed to determine 
whether any of the four boundaries (North, South, East 
and West) is a free-boundary. If this test is positive, 
the free-boundary conditions are incorporated. 

Namely, the program calculates here, for 
any free-boundary, the proper densities and velocities 
as dictated from the small-wave theory (for the 
lateral velocity components) and from the isentropic 
formula (for the axial velocity component). These formulas 
have been discussed in detail in Section 3. 

Next a check is made to determine if ISTEP=0 which, if 
positive, results in an iteration on the hydrodynamic and 
turbulence equations to improve the starting profile. 

After this series of tests, a call is made to STRID3 which 
yields the velocity and pressure fields. A call to STRID4 
performs similar calculations for all other dependent 
variables . 

Then the program calculates and prints out the local values 
of wall fluxes at every existing wall and their averages. 
Useful quantities calculated here include mean velocity and 
dynamic head, total areas, wall shear stresses, drag 
coefficients and Stanton numbers. 

Finally the value of ISTEP is incremented and, depending 
upon a pre-determined criterion, calculations are continued 
for the next station or terminated. This criterion is the 
limiting value of the forward distance. If ZU < ZLAST, 
control is returned to a point in MAIN just after the CALL 
BEGIN statement, and the process repeated for the next 
station. 

4.5 Detailed List of Program Variables 

All variables used in common statements are defined in 
Appendix d. Variables which appear locally are defined 
by their use in the computer program and are not defined 
in this report. 
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5. THE INVARIANT PORTION OF SHIP 

5.1 Introduction 

Many of the calculations do not need changing for different 
boundary conditions, etc. In this chapter the function of 
the invariant portion of SHIP is given. Again, Important 
variables are defined. 

5.2 STRIDA 


STRIDA consists of three separate sections which are called 
through ENTRY statements. STRIDO and STRIDl, the first two, 
calculate quantities related to the grid arrangement. The 
information supplied to STRIDA through BLOCK DATA is LCV, 

MCV and ZETA(I), ETA(J). The former represent the number of 
({)-control volumes along the ^ and n directions; and the latter 
denote the grid disposition in non-dimensional coordinates. 

Given the above information, STRIDO, the first member 
subroutine of STRIDA, computes the maximum niimber of grid 
nodes in the I and J directions. This is done in the 
following sequence: 


L 

= LCV 

+ 

1 

M 

= MCV 

+ 

1 

LPl 

= L + 

1 


MPl 

= M + 

1 



It is emphasized here that users must ensure that LPl and 
MPl are always less than or equal to IMAX and JMAX 
respectively. The latter represent the maximum dimensions 
of all variables in the respective directions and are given 
values accordingly in BLOCK DATA, Following the above 
sequence, the integer arrays JM and NFM are filled in 
accordance with equation (4.3-3). 


The second member of STRIDA is STRIDl. STRIDl is called 
once at the beginning from MAIN and never again if the 
calculation domain cross-section remains unchanged. For 
domains of axially-varying cross-sections it is, however, 
called at every step from MAIN. It is in STRIDl that the 
physical coordinates x and y are computed from the values 
of AGEOM, ETA and ZETA, as follows: 

X(I) = ZETA(I)*BXE (52 

Y(J) * ETA(J)*BYN 
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ETA(J) is initially calculated as: 

ETA(J) = ETA(J-l) + .AGE0M**(J-2.)*DELY (5.2-3) 

where : 

DELY = BYN* (1.0 - AGEOM)/(1.0 - AGEOM++M) (5.2-4) 

Here, BXE and BYN represent the width and height of the 
calculation domain respectively ({Xg - and {Yj^ - Yg} 

respectively); as such, they are to be specified oy the 
user in BLOCK DATA. It should be mentioned that any ETA 
distribution can be used; the distribution (5.2-3) is just 
an example; it makes use of a geometric series for specifying 
successive intervals between ETA's. Larger values of AGEOM 
result in the grid being "crowded" closer to the y = 0 
surface. 

The width and height of the calculation domain are allowed, 
in SHIP, to be varied with axial distance; their new values 
being calculated by the expressions: 

BXE = BXEU + DBXEDZ*D74 

BYN = BYNU + DBYNDZ*DZ (5.2-5) 

where BXEU, BYNU represent the upstream (i.e. of previous 
station) values of BXE, BYN and the quantities DBXEDZ, DBYNDZ 
(i.e. the slopes of BXE and BYN respectively) are specified 
by the user. 

The ratio of cell areas in the ^-p plane between upstream 
and downstream stations is also calculated in STRIDl, and 
stored as ARAT; the expression for this ratio is: 

ARAT = (BYNU*BXEU)/(BYN*BXE) (5.2-6) 

The subsequent operations in STRIDl deal with calculation of 
various inter-grid distances required later in the 
calculation of the coefficients. These distances are 
shown in Figure (11) , which illustrates the grid in the 
n-5 plane and the nomenclature described below. 
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The dashed lines in figure' (11) .join the control-volume 
faces normal to the x and y directions. These are the 
control volumes for all dependent variables except the u 
and V velocity components. The con trol-volume faces pass 
mid-way between the grid nodes except near the boundaries 
where they pass through the boundary grid nodes. Thus, 
the control-volume faces always pass through points where 
the velocity component normal to the faces are stored. A 
normal velocity component at a control-volume face is 
presumed to prevail over that whole face. 

The control-volumes for each of the velocity components 
u and v are displaced along the directions of these 
velocities. The control-volume faces normal to each of 
these directions pass through grid nodes on either side of 
the velocity component in question. Figure (12) illustrate: 
this point. 



V{ I + l, J) 


FIGURE 12: 


U- AND V- VELOCITY CONTROL VOLUMES 



The independent variables x, y are stored as the coordinates 
X(I) and Y(_J) of the grid lines. Certain related quantities 
are necessary and are calculated and stored here, either 
once-for-all if the domain cross-section remains unchanged, or 
at every station if it changes. Among these XS(I),. XDIF(I) 
and XSU(I) are derived from the values of X(I) such that: 

XS(I) denotes the x-direction length of a main control- 
volume around the node (I,J);XSU(I) is the x-direction 
length of the control-volume for U(I,J). (Note that XSU(I) 
is the same as XDIF(I) except near the boundaries of the 
calculation domain). Incidentally XS(I) is also the 
distance between the locations of U(I,J) and U(I + 1,J) and 
performs the work of XDIF(I) for the equation for u. The 
meanings of YS(J), YDIF(J), and YSV(J) should now be 
obvious. The geometric quantities associated with the grid 
system are defined in Table 5. 

Because the faces of the main control-volume are defined 
to pass midway between the grid nodes, interpolation factors 
for the calculation of variables on these faces are equal 
to 0.5, except near boundaries where they are either 0 or 1. 
For the control-volumes appropriate to U- and V- velocity 
components however, the interpolation factors can differ 
from 0.5 if the spacing between grid nodes is chosen to be 
non-uniform. For this reason, interpolation factors are 
calculated in STRIDl and stored as FXP(I), FXM(I), FYP(J) 
and FYM(J) , The subscript refers to a grid node. The value 
of U for example at a grid node (I,J) is given by: 

FXP(I) * U(I + 1,J) + FXM(I) * U(I,J) (5.2-7) 

It is obvious from the above that FXM(I) is simply 

1.0 - FXP(I), The expressions for FYP and FYM are similar. 
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TABLE 5. Definition of Grid Geometric Quantities 


NO 

GRID 

' >QUANTITy 

MEANING 

• ; 1 

1 . 

. I 

' X(I) . 

Physical coordinate in the ^-direction; 
represents x in cartesian coordinates. 

2 

.. XDIF(I) 

.. i.i 

The difference between X(I) and X(I-l); 
used as the distances 6 in calculating 
^-direction diffusion flux of 4): 
r^A(j)/6* 

3 : 

XS(I) 

The ^-direction width of a main control 
volume. 

4 . 

XSII(I) 

The ^-direction width of a U-velocity 
control volume. 

5 

y(j) 

Physical coordinate in the n-direction; 
•represents y in cartesian coordinates. 

6 

YDIF(J).'I 

The difference between Y(J) and Y(J-l); 
used as the distances in calculating 
n-direction diffusion flux of 4>. 

7 

YS(J) 

i The n-direction width of a main 
control volume. . 

8 

YSV(J) 

The n-direction width of a v'-velocity 
control volume. 



i 
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Calculation of the quantities tabulated above completes 
the tasks performed by STRIDl. STRIDl is called once 
from MAIN to calculate the Initial grid quantities; 
thereafter, if lARCH =2, i.e. variable areas of the 
integration plane are to be accounted for, STRIDl is called 
from MAIN at every s^ep, after a call to GEOMOD has 
provided the necessary information for the cross-section 
change. 

The third section of STRIDA, namely STRID2, aids partly in 
the calculation of the coefficients in the finite-difference 
equations. STRID2 is used to calculate and store the 
convective mass velocities (i.e. the terms in square 
brackets in equation ( 2 . 3 - t 2 ) crossing the control volume 
faces along the £- and n-directions. 

The arrays GX and GY are used, respectively to store these 
values. STRID2 is called for each Integration plane, once 
at the beginning of STRID3, then again from STRID3 after 
the grid has been adjusted to fit the downstream geometric 
configuration. Finallv STRID2 is called from STRID4, after 
the u and v velocities have been corrected. 

5.3 STRIDE 


STRIDE is the largest and perhaps the most important 
subprogram. STRIDE is the main machinery of the SHIP 
program and contains instructions to calculate the 'A' 
coefficients in the finite-difference equations (3.2-4) 
and (3.2-7) to (3.2-9). The coefficients are assembled 
in STRIDE based on the convection and diffusion fluxes 
in the three coordinate directions. STRIDE contains 
two member subroutines: STRIDE and STRID4. 

It is in STRIDE, that the finite-difference equations 
appropriate to U, V, W and PP are assembled, in that 
order; while STRID4 assembles the coefficients for all 
other equations such as for kinetic energy of turbulence, k, 
etc. The first entr 3 »^ statement to STRIDE is RTP.ID3 . 

STRIDE calls first STRID2 to assemble the mass fluxes GX, 

GY and GZ crossing the respective cell areas. These values 
are computed at any section, from the upstream distributions 
of u, V and w. 

The first equation assembled by STRID3 is the u-inomentum 
equation. Later, in sequence, the v, w and pressure- 
correction equations are solved. Since the structure of 
FORTRAN instructions is similar for all equations, it is 
sufficient here to explain that for one of them only. 
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For the u-momentum equation, the details are as follows: 

The statements defining ISTB and JSTR specify the 
starting I value and J value for the equation. For the 
u-velocity equation, ISTR is 3 and JSTR is 2; for the 
v-velocity, ISTR = 2 and JSTR = 3; for others 

ISTR = JSTR = 2. The next imnortant instruction 
is a call to GAMMA; GAMMA calculates the diffusion 
coefficients that enter the finite-difference equations. 
GAMMA in turn, makes a call to GAMOD whereby account is 
taken of the modifications for the -boundaries. In the 
present version of SHIP, the call to GAMOD is redundant. 

A subsequent call to SOURCE computes the source terms in 
the equations. It should be noted that the SU and SP are 
such that : 

S » SU +SP*(j)p 

where S is the source term of the equations. This source 
term is then integrated in STRID3 over the control 
volume. 

The calculation of the A coefficients starts first with the 
calculation of the upstream convection coefficients, denoted 
by ALZ(IJ). These are merely the mass fluxes over the cell 
areas normal to the z direction. In the cross-stream the 
coefficients are calculated first for the near-boundary 
nodes. 

The first statements calculate coefficients arising from 
the J = 1 boundary (i.e. pertaining to AYM for J = 2) ; 
and the next statements calculate those for the 1=1 
boundary. Then in a DO loop, the coefficients for the 
other interior nodes are calculated. The AL-s represent 
the convection part (the L*s in equation (3.2-5)), and the 
T-s represent the diffusion part (the hyphen stands for 
X or Y depending on the direction considered) . The high- 
lateral-flux modification is then introduced. Then, the 
pressure gradient term is calculated; DU is the area over 
which the pressure forces act; and SU is the source term 
summed up along with the contribution of others. When all 
the coefficients are assembled, a call to SOMOD is made to 
determine whether any of the earlier computed coefficients 
and source terms' require modification, to account for the 
particular flow-geometry. Once this is achieved, a call 
to SOLVE permits solution of the linear algebraic 
equations to be obtained by sweeps of TDMA. 

The other subsections of STRID3 perform identical operations. 
First the v and w equations are solved in a manner similar 
to that for u. After this is done, the requirements of 



overall mass continuity are examined. Then, the coefficients 
for the pressure-correction equation are assembled. The 
calculation of these coefficients is similar to that for 
other equations; but it is to be noted' tha^ for the pressure- 
correction equation, AZ is identically zero; and the term SIJ 
is the mass source. It is important to note here that- 
the mass flux from the downstream face (i.e. pw+) is 
calculated from an updated density at the downstream face. 

The pressures and velocities are subsequently corrected 
and this concludes STRID3. 

The assembly and solution of the finite-difference 
equations of other dependent variables is done in STRID4. 
Before STRID4 is started, a call is made to f3TRID2; this is to 
calculate the new values of GX, GY based on the newly 
calculated u and v velocities. The operations in STRID4 
are similar to those explained in STRID3. 

At the end of STRIDE, the directions to perform the TDMA 
sweeps are given. The indices IXY, ISWP and JSWP are 
changed through the statements 

IXY=2-IXY 

ISWP=3-ISWP (5.3-1) 

JSWP=3-JSWP 

The completion of STRID3 and STRID4 yields new distributions 
of all variables and the pressure p at a downstream 
location. The repetition of STRID3 and STRID4 at 
several forward steps is controlled by the main program 
SHIP. 


5.4 Solution Procedure for Algebraic Equations 
The Subprogram SOLVE 

The function of the subprogram SOLVE is to arrange for the 
solution of the finite-difference equations, for each 
dependent variable NV, to be obtained. The solution 
procedure used is the application of a pair of TDMA 
traverses, one in each of the C-and n-directions . 

SOLVE has three major sub-divisions. The first ends with 
the statement: 10 CONTINUE. It is in this part that the 
finite-difference coefficients are assembled in readiness 
for the subsequent operations. It is here that our fully- 
implicit difference scheme is implemented during the 
coefficient-assembly process; incorporation of other 
schemes will require some modification of this part of SOLVE. 
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The second part of SOLVE is concerned with TDMA traverses 
in the ^-direction and ends with statement: 21 CONTINUE 

The third and final part of SOLVE starts with this statement 
and concerns the TDMA traverses in the p-direction. 

A call to SOLVE is made from STRID3 and STRID4, once for 
each dependent variable NV. This call achieves TDMA 
traverses in both the and p-directions ; however, 
which traverse is made first depends upon the value of the 
index IXY. The first traverse direction will be E, if 
IXY = 1, and p if IXY = 2. IXY is set alternately to 1 or 
2 by the statement: IXY = 3 - IXY which concludes both 
STRID3 and STRID4. 

Along each traverse direction, the direction of sweep i.e. 
whether from the first grid node to the last or vice versa, 
depends upon the value of an index, ISWP in the ^-direction 
and JSWP in the p-direction. Each of these indices takes 
on a value 1 or 2 by statements which follow the end of 
parts 2 and 3 of SOLVE. For example, the statement 
following: 21 CONTINUE, reads JSWP = 3 - JSWP. A value 
1 implies that the sweep direction is from first to last 
grid node and 2 implies vice versa. 

On occasion, it may be required to perform more than one 
pair of TDMA traverses in a ? - p plane, for any given 
dependent variable. The number of pairs of TDMA traverses 
is set by values assigned to an index array NSWP(NV). The 
program is set up with NSWP values equal to unity except 
for the pressure-correction equation for which NSWP(NPP) = 3. 
For locations near injection points the NSWP values are 
increased to ensure stability and accuracy. 

5.5 ■ Printout of Field Values of Dependent Variables 

PRINTCISKIP. JSKIP) 


It is frequently required to print out the contents of the 
F- array. This task is performed by subroutine PRINT. 

PRINTCISKIP, JSKIP) provides for the printout of F(I,J,NV) 
where NV can attain a maximum value of NFPMAX, which is set 
as : 


NFPMAX = NNV + 3 (5.5-1) 

The three extra values representing NVP, NRO and NMU, i.e. 
pressure, density and effective viscosity respectively. 

The decision as to whether a particular variable NV is 
printed out or not, depends upon whether the corresponding 
IPRINT(NV) is equal to unity or not. 
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The printout of each dependent variable NV is given a 
heading stored in TITLE( . . . ,NV) . The formal parameters 
ISKIP and JSKIP permit the selective skipping of columns 
(I) and rows (J), when it is not required, for any reason 
to printout the complete array of values of each variable 



6. THE PHYSICAL MODELLING SECTION OP SHIP 

6 . 1 Introduction 

This chapter describes the portion of SHIP which calculates 
the physical properties required in the solution procedure. 

6.2 AUX 

The general policy is to confine all tasks associated with 
physical modelling to subprogram AUX. Thus the calculation 
of density, p, effective diffusion coefficient r, and sources 
and sink's S of the dependent variables are performed in 
separate member subroutines in AUX, namely, DENSTY, GAMMA, 
SOURCE, VISCOS and SPECIE. 

AUX to a large extent, is an invariant subroutine; it need 
not be changed unless different physical laws and turbulence 
model need to be incorporated. AUX is not called as such 
from any subroutine; instead the various parts mentioned 
above are referenced to when necessary. 

6.2-1 DENSTY 


This subroutine calculates the densities pp y {RHOCI,J)} 
and Pp {RHODCIJ)} in accordance with the discussion of 
Section 3.2-6. At each forward step new RHOD(I,J) are 
calculated from guessed values, of pressure. After 
calculating the starred velocity fields, RHOD(I,J) are 
changed to correspond to the new velocity distribution. 

6. '2-2 VISCOS 


The function of VISCOS is to calculate the viscosity, 
both laminar (if necessary) and turbulent. In its present 
form, the laminar contribution has been neglected, since 
the flows under consideration are highly turbulent. Hence, 
the major function of VISCOS is to store the turbulent 
viscosities . in the array AMUT(I,J). 

VISCOS also performs the function of calculating for the 
N,S,E and W boundaries, the effective boundary diffusion 
based upon the semi-logarithmic law-of-the-wall . These 
diffusion coefficients are stored respectively in the 
arrays GAMN(I), GAMS(I), GAME(J), and GAMW(J), which are 
later substituted in appropriate GAM locations by GAMOD. 
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6.2-3 GAMMA 

This part of AUX is used to set values to the exchange 
coefficients array GAM(I,J). These are calculated simply 
as: 


GAM(I,J) = AMUT(IJ)/PR(NV) (6.2-1) 

A call to GAMOD is then made in order to permit any 
modifications to be made to the GAM array. 

6.2-4 SOURCE 


The/ section headed by SOURCE is concerned with the 
calculation of the finite-difference form of the source 
terms in the equations stated in Table 1. The 
source terms for each variable are programmed under 
separate subsections; and choice is made through the NV 
indices (e.g. NVU, NW etc,). The source terms calculated 
in SOURCE are for unit volume and represent volume averages; 
they are multiplied in STRIDE by the volume of the appropriate 
control volume. The production term for the kinetic-energy 
of turbulence, k is calculated also in AUX, and is stored 
in the array GENR(IJ). 

It may be noted that all source terms contain two components 
SU and SP. SU is that part of source terms which is 
calculated completely from upstream values of the variable; 

SP is the linearised part. in total they express combinedly 
the relatioh 

S^ = SU(I,J) + SP(I,J)*F(I, J,NV) (6.2-2) 

It is necessary to ensure that SP is always negative . 

6.2-5 SPECIE 


Subroutine SPECIE is used to calculate the species mass 
fraction from knowledge of the element mass fractions, 
pressure and temperature. Appendix A gives a full 
development of the technique employed. 
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• PROBLEM^DEPENPENT SECTIONS OF SHIP 


7.1 Introduction 

This chapter completes the detailed description' of the 
functions of the various portions of SHIP by discussing 
the subroutines which must be changed in order to specify 
a particular problem. 

7.2 General Policy 

As mentioned briefly earlier, the main machinery of the 
numerical calculation procedure is devoid of any problem- 
specification information; the latter is provided through 
subroutines BLOCK DATA and ALLMOD. It is envisaged 
therefore, that detailed modifications to the above- 
mentioned routines, may require to be made for each new 
problem and should be done so with care. 

7.3 BLOCK DATA 

BLOCK DATA serves to provide values to fluid properties, 
grid distributions, program control parameters and other 
information specific to each new problem, via DATA 
statements. The use of BLOCK DATA permits the program to 
be run with compilers common to both CDC and IBM machines. 

In the present version of SHIP the above information is 
given in chapters. Thus, Chapter 1 deals with preliminary 
information such as SMALL, GREAT etc. Chapter 2 provides 
the information necessary to specify the grid and geometry. 
Chapter 3 is concerned with the dependent variables 
information and Chapter 4 with the physical properties data. 
Chapter 5 provides some starting values whilst Chapter 6 
is concerned with step control. Chapter 7 provides the 
fixed boundary conditions and indices; finally Chapter 8 
provides the indices required for printing-out. 

7.4 ALLMOD 


Subroutine ALLMOD contains instructions for incorporating 
specific information regarding the flow geometry, and for 
any specific changes to the coefficients in the finite- 
difference equations, or to the variables themselves. 
ALLMOD is divided into several " MOD if ication" routines 
which will be described below. More specifically, it is 
composed of seven member subroutines BEGIN, GAMOD, GEOMOD, 
SOMOD, UPSTRM, INJMOD and INJMOT. 
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7.4-1 BEGIN 


The primary purpose of this subroutine is to provide 
initial values to all the dependent variable arrays, 
fluid— property arrays and other auxiliary arrays, so 
that the marching integration can start. The inlet values 
of velocities, temperatures and turbulence quantities are 
specified here. Also here, the turbulent Prandtl 

number for dissipation of kinetic-energy of turbulence 
is calculated from the relation: 


K‘ 


(C2-Ci)Cj 


(7.4-1) 


Experience with the k-e model of turbulence has shown that the 
above relation, though it applies strictly to the wall region, 
produces satisfactory results for regions far from the wall. 
This relation is derived by eliminating the convection terms 
from the conservation equation of e and by evaluating the 
production terms of e from the logarithmic law of the wall. 


The secondary purposes include provision of 'fixed' boundary 
conditions on the four boundaries of the calculation domain, 
and calculation of some auxiliary information required in 
the initialising process. 


The preliminary calculations for incorporating the free- 
boundary conditions are also made here, depending on the 
values of the indices KBCE, KBCW, KBCN, KBCS. A free boundary 
is implied when the value of the above coefficients is 3. 


It should be noted that the quantity MW/RT used in the above 
calculations for the free-stream, assumed MW = 28.93. This 
should change if the free-stream is not pure air. The user 
will be reauired to give careful attention to this section. 

7.4-2 GAMOD 


GAMOD is provided to allow the exchange. coefficients GAM(I,J) 
computed in AUX to be modified, as necessary, to incorporate 
the required boundary conditions. 

As a matter of fact, the contents of. GAMOD in the present 
version of SHIP are redundant. The required modifications 
are made directly to the 'A' coefficients in SOMOD. 

7.4-3 GEOMOD 


GEOMOD is the subroutine where modifications to the geometric 
configuration of the calculation domain can be prescribed. 


63 


To activate GEOMOD, I ARCH is set equal to two; otherwise, the 
geometry will be a rectangular solid set by the initial condi- 
tions. The code can handle any geometry in which there are no 
sudden or discontinuous changes in domain size and the cross 
section of the domain normal to the main flow direction is 
rectangular. The geometry is determined by four constants 
whose magnitudes are set in GEOMOD. These constants (with 
definition) are: 

DXWDZ = rate of change of distance between reference plane and 
west boundary with respect to distance in the main flow 
direction . 

DYSDZ = rate of change of distance between reference plane and 
south boundary with respect to distance in the main 
flow direction. 

DBXEDZ = rate of change of width of domain with respect to 
distance in the main flow direction. 

DBYNDZ = rate of change of height of domain with respect to 
distance in the main flow direction 

The values of these four quantities are determined by the given 
domain and can either be constants or functions of distance in 
the m.ain flow direction. 

7.4-4 SOMOD 


The purpose of SOMOD is to facilitate modifications to the 
source terms SU(I,J) and SP(I,J) and the coefficients in 
the finite-difference equations, AXM(I,J), AXP(I,J) etc. 
SOMOD, together with GAMOD, achieves the final form of the 
finite-difference equations before they are solved. At 
present, the function of SOMOD is to incorporate the wall 
functions and the symmetry boundary conditions; but a 
greater amount of modification, such as creating internal 
obstacles etc. , can be achieved through SOMOD. SOMOD is 
divided into several sections each dealing with a separate 
variable. For the equation for pressure-corrections no 
modifications are made. 

As mentioned earlier, the provision of wall-functions is made 
through SOMOD. This is achieved as follows: first, 

coefficients linking boundary grid nodes with their immediate 
neighbours inside the calculation domain are set to zero for 
each variable. This is done because, by incorporating the 
wall functions, we are prescribing the values of the fluxes 
directly; and hence whatever values are calculated earlier 
in STRIDE for the coefficients must be set to zero. Then if 
an index, for example, KBCS for the ^outh boundary, is set 
equal to unity, the appropriate wall flux is calculated and 
fed in through SU and SP. This calculation uses the value 
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of GAMS(I) (calculated in AUX, and brought here through a 
COMMON statement); and the corresponding flux, FLUXS(I,N1^) 
of the dependent variable is stored for purposes of print- 
out. If the index is other than unity, no change is made 
to the coefficients. Similar indices KBCN, KBCE and KBCW 
are used for the North, East and West walls. Note that 
for the v-velocity, for example, there are no modifications 
for the North and ^outh walls since they are noi'mal to that 
velocity. For the w-velocity the wall functions for the 
above wall boundaries are similar to that for u. Such 
checks and modifications are made for the variables u, v, 
w, k, h and f. 

The modifications to the source term for k consist of 
altering the generation and dissipation based on shear-stress 
from the wall functions. For the dissipation equation, 
the value of dissipation is fixed according to equation 
(3,4-4), by modifying SU and SP. Setting SP to a large 

30 

negative value (-10' ) and SU to a large positive value 

multiplied by the required value to be fixed ), 

we nullify the effect of other coefficients (since 
they will be divided by SP) , and obtain e equal to the 
desired value This practice is adopted in general 

for fixing any variable at a desired value. 

The user of course may be required to provide other types 
of boundary conditions himself, again through SU and SP. 

In this respect, choices of the index values other than 
unity can be profitably put to use. 

7.4-5 UPSTRM 


The subroutine UPSTRM makes provision to store upstream 
(i.e. previous integration plane) values of pressure which 
are necessary in calculating the source term of the axial 
velocity component. The upstream values of the width 
and height of the integration domain are also stored here, 
for variable geometrical configurations, UPSTRM is also 
a convenient place to calculate the axial direction mass 
flow rate. It should be noted that upstream values of any 
other variables apart from pressure are not required to 
be stored. This should be the case only if iterations 
were necessary. 

7.4-6 INJMOD 


This subroutine is primarily used to specify the correct 
flux of a given variable through the array FINJ(NV). This 
array is used to modify the source term in SOMOD when 
injection occurs through the ^outh wall. 
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7.4-7 INJMOT 


This subroutine has a function identical to that of INJMOD 
above, but for injection occuring through the Noi'th wall. 
INJMOT and INJMOD may be called either simultaneously or 
separately from MAIN; this depends on whether injection 
occurs through both walls or through either wall and is 
also decided in MAIN. 
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8. SOME USER’S GUIDANCE IN ADAPTING SHIP FOR A GIVEN PROBLEM 

8.1 General Remarks 

The following steps are necessary in adapting SHIP for a 
given problem. 

a. Provide the following information in BLOCK DATA 

• Program control parameters 

• Grid and geometric specifications 

• Printout control parameters 

• Solution procedure parameters 

• Fluid property values 

• Thermodynamic and hydrodynamic data appropriate 
to the problem 

• Fine-tuning of the turbulence model constants 
to optimise the results. 

b. Provide adequate initial values to variables in BEGIN. 

It may be required, for example, to chose a grid 
disposition (ZETA (I) and ETA (J)) to suit known initial 
distributions of dependent variables. 

c. Check if known boundary conditions are correctly set 
and unknown values correctly calculated. 

d. Incorporate the exact geometric configurations of the 

problem into the appropriate sections of ALLMOD. 

e. Arrange in PRINT, for requisite printout 

8.2 List of Input Variables 

Following is a list of input variables to be supplied to 
the SHIP program; they correspond to the variables defined 
in the BLOCK DATA routine at the beginning of the program. 
The following information is given for each variable. 

(i) the FORTRAN symbol 

(ii) the meaning 

(iii) for some variables, the recommended value 

is given in the last column. 
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SYMBOL 



Por Program Control I 


They correspond to the dimensions of the 
major arrays. 

The number of variables (NV) for which the 
F array is dimensioned. 

They are variable numbers having values 
ranging from 1 to 10 . 

The step number; it is initialised in BLOCK DATA. 0 

Specifies whether the flow is to be treated as: 
one of constant geometry, ^ 

or varying with axial direction. 2 

The initial step length. 

The maximum step length. 

Expansion factor for step length. FRA at step n, 

(FRA ) , is equal to FRA .*EX or FRAM whichever 
is smaller. 

Small and large numbers used at various places 10 

in the program, 10 

2. Grid and Geometry Variables 

Total number of control volumes in the ? 
direction. 

Total number of control volumes in the n 
direction. j 

The disposition of grid nodes in the ? 
direction. ZETA ranges from 0 to 1 always 
and is given, in general, by 

zetaci) . Sii : 


The locations of the grid nodes in p direction 
ETA(J) 
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SYMBOL 

MEANING - 

REC 



VALUE 


KBCN 

Specify the nature of the North, South, East 


KBCS 

and West Boundaries# Thus for: 


KBCE 

Wall boundary 

1 

KBCW 

Symmetry axis 

2 


Free boundary 

3 

BYN 

Height of the domain. BYN must be 
initialised to the value at the inlet plane. 


BXE 

The extent of solution domain in the 
x-direction. 


ZU, ZD 

Axial locations of upstream and downstream 
computational planes - must be initialised 



here. 

0.0 

DZ 

Axial step length. The value set here is 
overwritten by that calculated from FRA 
specification . 


DYSDZ 

dy /dz - initialised here to zero. 

O 

0.0 

DBYNDZ 

d(BYN)/dz - initialised here to zero. 

0.0 

ZLAST 

The z location at which solution is terminated. 
1 3. Control Variables for Printout 


ZRE 

The 'z' locations of the axial stations at 
which the distributions of the variables 
are printed out. 


LASTED 

The total number of axial steps at which the 
variables are to be calculated. At present it 
should be set to a large value so that the 
solution is terminated by ZLAST. 


NPJUMP 

The number of steps between successive 
printouts of the variables. 


I C JUMP 

Has the same meaning as NPJUMP; but it controls 
printout of wall fluxes. 


IPRINT 

Specifies whether the distribution of the NV 


(NV) 

variable should be : 



Printed 

1 


Not printed 

0 








4. Solution Procedure Parameters 

The value of ISOLVE(NV), specifies whether 
or not variable number NV is to be solved 
for. 

ISOLVE = 0 does not solve for the variable 
ISOLVE = 1 solves for the variable. 


Number of TDMA sweeps performed on variable 
NV at each step. 


6 for NV=1 
and 3 for 
all others 


These variables control the order of 
TDMA sweeps at each step. Recommended 
starting values are unity for all. 

5. Turbulence-model Constants 


1 , 1 . 1 . 


PRLAM 

(NV) 


Laminar Prandtl number for variable NV 


1.44 

1.92 


0.09 


K in the log law. 

E in the log law. 

Turbulent Prandtl numbers for the variables. 
Their values are set initially to unity. 
However for the dissipation of e, PR is later 
calculated from the expression 
PR(NV) = k2/((C2-Ci) Cpi); further, PR for 
enthalpy is somewhat less than unity, equal 
approximately to 0.90. 

Factor relating turbulence energy to mean 
motion energy. 

The P function in the "wall functions” 
(resistance of laminar sublayer). 

6. Fluid Properties 


0.42 


9.00 


0.004 
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SYMBOL 

MEANING 

[H33HI 

nRSMI 

WM 

GASCON 

Molecular weight and universal gas constant 
for the gases considered, At present GASCON 
is taken to be 8314 m^/sec^OK 


AMUREF 

Reference laminar viscosity. 


DEN 

Reference density. 


CP 

Specific heat. 


CPDCV 

Specific heat at constant pressure divided by 
specific heat at constant volume. 

1.4 

HO 

Enthalpy of formation. 



7. Starting Values 


a 

Inlet X- , y- and z- direction velocities. 


PIN 

Inlet pressure. 


TIN 

Inlet free stream temperature. 


RETRAN 

Transition Reynolds number. . 



8. Boundary Conditions 


IINJ 

IINJT 

Indices to control entry to INJMOD and INJMOT, 
respectively . 


INJSTP 

INJTOP 

Values of ISTEP at injection locations on 
the South and North walls, respectively. 


ZINJ 

ZINJT 

z locations of injection, for the South 
and North walls, respectively. 


TINJ 

TINJT 

Jet temperatures, for the South and North 
walls, respectively. 


FLOINJ 
FLO I NT 

Mass flow rates injected by a jet through 
the South and North walls, respectively. 


I ADI AN 
^ lADIAS 
I ADI AW 
lADIAE 

k 1 

Indices specifying whether the North, 
South, West and East walls, respectively 
are adiabatic (=1) or under constant 
temperature (=0). 
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SYMBOL 

MEANING 

REG 

. . VALUE . . 

TWALLN 

TWALLS 

TWALLW 

TWALLE 

North, South, West and East wall 
temperature, respectively. 
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The units used by the program are: Kg, m, sec. It may be 
convenient, however, for a number of input quantities to be 
specified in other units; this can be done provided that 
the conversion to Kg, m, sec is then effected internally. 

A comnlete list of PORTRAN variables is given in 
Appendix C. 

8.3 The Grid 

• The axial grid is specified via the variables FRA, FRAM 
and EX, defined in the list above. 

• To change the grid distribution in x and y, while 
retaining the same total number of nodes in each direction, 
simply reset the ZETA (I) and ETA (J) arrays. 

• If the number of grid nodes is to be changed, then LCV 
and MCV must be reset, and ZETA and ETA must be 
respecified. Note that the number of ZETA values provided 
must be LCV + 2, and the number of ETA values must be 

MCV + 2. 

• When the number of grid nodes is increased it may be 
necessary to redimension some of the arrays, as it is 
described below. 

• Note that, whenever the grid distribution is changed, 
it will be necessary to reset the inlet conditions 
accordingly. 

8.4 Dimension Changes 


The quantities IMAX, JMAX and NNV must be set to correspond 

to the reauired dimensions of the arrays in I, J and NV 
(thus IMAX > LCV + 2, and JMAX > MCV + 2). The arrays must 
then be dimensioned as given in'*the table below, and 
the EQUIVALENCE statement must be modified to read: 
EQUIVALENCE ( F( 1) ,PP( 1) ) , (F(IMAX*JMAX + 1), U(l)), 

(F(2+ IMAX* JMAX + 1), V(l)), ( F( 3* IMAX* JMAX + 1), 

W(D), (F( 4* IMAX* JMAX + 1), AKE(l)), (P( 5* IMAX* JMAX + 1), 
ADS(D) . 
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Array 


Must be dimensioned 


1 . General 
F 

P, EHO, GAM, SU, SP, DU 
DV, DW, AXP, AXM, AYP, 

AYM, A2, GX, GY, GZ 

CY 

A, B 
CYU 

DRHODP, PU, RHOD 

A 

PP, U, V, W, AKE, ADS 
NFM 

ZETA, X, XS, XDIF, XSU 
YSR, ETA, Y, YS, YDIF , YSV 
FXP, FXM 
FYP, FYM 
JM 

I SOLVE 

RELAX 

NSWP 

IPRINT 

TITLE 

ZRE 

AMUT, CPMN 

PRLAM, PR 
FLUXN, FLUXS 

FLUXW, FLUXE 


( IMAX* JMAX*NNV) 
(IMAX+JMAX) 

At least (LCV + 2) 

At least the greater of 
(LCV + 2) and (MCV + 2) 

At least (LCV + 2) 

(IMAX+JMAX) 

( IMAX* JMAX) 

At least (NNV + 7) 

At least (LCV + 2) 

At least (MCV +2) 

At least (LCV + 2) 

At least (MCV + 2) 

At least (MCV + 2) 

At least (NNV) 

At least (NNV + 3) 

At least (NNV) 

At least (NNV + 7) 

(6, at least NNV + 7) 

At least (no, of printout 
stations required) 

At least (IMAX* JMAX) 

At least (NNV) 

At least (LCV + 2, NNV) 

At least (MCV + 2, NNV) 



Array 


Must be dimensioned 


GAMS, GAMN, 


At 

least 

(LCV 

+ 

2) 

GAME^ GAMW 


At 

least 

(MCV 

+ 

2) 

2. Subroutine 
XI 

PRINT 

At 

least 

(LCV 

+ 

2) 

X2 


At 

least 

(MCV 

+ 

2) 

3. Subroutine 

AUX 






DUIDXJ 


(3 

,3) 




GENR 


At 

least 

(LCV 

+ 

2, 

4. Subroutine 

ALLMOD 






FINJ, FINJT 


At 

least 

LCV 




8.5 Interpretation of Output 
( a) Preliminary output 

This is provided by the program before the marching 
integration commences. It comprises the conditions at 
inlet, printed out as described below. 


The variables printed are as follows:* 

• the x-direction velocity component (m/sec) 

• the y-direction velocity component (m/sec) 

• the z-dir.ection velocity component (m/sec) 

2 2 

• the turbulent kinetic energy (m /sec ) 

• the dissipation rate of turbulence energy (m /sec ) 

• the stagnation enthalpy (Kcal/Kg) 

• the hydrogen mass fraction in any form 

• the thermodynamic temperature (°K) 

2 

• the static pressure (N/m ) 


♦Note that the user is free to suppress printout of 
any of these variables by use of the IPRINT array. 
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the 


density (Kg/m^) 

• the diffusion coefficient for ^(Kg/iri sec) 

• the mass fraction of Hg 

• the mass fraction of Og 

• the mass fraction of OH 

• the mass fraction of HgO 

• the mass fraction of H 

• the mass fraction of 0 

• the mass fraction of Ng 

• the mass sources in the field 


The printout is of the same format in each case. First is 
printed the heading (i.e. TITLE (6,NV)) for the variable in 
question. Then the complete two dimensional field is 
printed, with I arranged horizontally, and J vertically. 
Also shown are, at the far right hand side, the Y-direction 
grid locations for each J, normalised by the height of 
the domain, and below the array, the X-direction grid 
locations for each I, in (m). 

(b) Main Output 

The output provided during the marching integration is 
of four kinds, as follows: 

(i) At every step, a single line of information 

is printed out. Variables are provided, which. have 
the following meanings: 

ISTEP - step number 

ENTRN - the entrainment in the Y-direction 
through the North boundary 

TOTM - the total mass flux in the axial 
direction 

TOTW - the total momentum flux in the axial 
direction 

TOTH - the total convective enthalpy flux 



T0TN2 

TOTF 

FLOBOT 

FTOP 


the total convective flux of Nitrogen 

the total convective flux of Hydrogen 
in any. form 

the jet mass flux through the South wall 

the total Hydrogen in any form mass flux 
at the North boundary 


The headings provided for the above variables are: 
ISTEP, E, M, MO, H, N2, F, B and BT respectively. 


(ii) At every step another single line of information 

is printed out. Variables are provided which have 
the following meanings: 

ISTEP - step number 

ZU - axial location, ^ (m) ’ 

Y(MPl) - the location of the north boundary 

(useful in checking the axial variation 
of the domain's height). 

SUM - the sum of mass sources, indicating 
the imbalances in the field. 

The headings provided for the above variables are: 
ISTEP, ZU, DEL and MASS IMBALANCE. 


(iii) Every ICJUMP steps, axially directed wall shear 
stresses are printed out; drag coefficients and 
Stanton numbers are also printed out. What is 
printed out is, for I from 2 to LCV + 1. 

• AREAS, AREAN, AREAW and AREAE - the areas 
(in m2) of the South, North, West and East 
walls adjoining the near-wall control volumes 
for the .1 in question 

• TAUS, TAUN, TAUW and TAUE - the local wall shear 
stresses at the South, North, West and East walls. 

• CFS, CFN, CFW and CFE - the local wall shear 
stresses at the South, North, West and East 
walls, divided by the mean dynamic head (ipw^) 
at the station under consideration 

• STANS, STANN, STANW,and STANE - the Stanton numbers 
at the South, North, West and East walls 

Apart from the above local quantities, are also 
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printed out mean quantities as follows: 


• WBAR - the mean value of w across each 

plane of calculation 

• DHEDTD - the mean dynamic head at each plane 

of calculation 


• CFSAV , 
CFNAV , 
CFWAV , 
CFEAV 


the mean values of CFS, CFN CFW 
and CFE defined above, respectively. 


• STASAV, 

STANAV, - the mean values of STANS, ST ANN, STANW 
STAWAV, and STANE defined above^ respectively, 

STAEAV 


It is also in this context that the following 
quantities are printed out for testing purposes: 


• FLOWT - the mass flow rate across each plane of 

calculation 


• DHED - the total momentum across each plane 

of calculation. 

(iv) Every NPJUMP steps, and/or as specified via ZRE, 

complete variable printout is provided as described 
in Section 8.5 (a) above. 



9. 


RESULTS AND DISCUSSION 


9 .1 Introduction 

Sample results for three typical cases of interest are presented 
in this cliapter. It should be noted that only the third case is 
the novelty of this section which has, otherwise, been taken 
from a previous report {6} of the HISS program. This case concerns 
parallel hydrogen injection in a domain whose North boundary wall 
expands with axial distance. The results for the other two cases 
have been taken from {6} and have been obtained by using HISS rather 
than SHIP. However, test runs for identical cases, using the latter 
program' proved that the results thus obtained are very similar to 
the ones reported here. 

9 . 2 Presentation of Typical Results 

9.2- 1 Introduction 

Results for ten cases were calculated for the contract {6} 
by the technique described in the first eight chapters. 

For each case profiles of 19 Variables described earlier 
were computed for a grid of 12 grid lines transverse to the 
flow and 20 grid lines normal to the flow for a total of 
240 points at each station in the main flow direction. 

Computations were made for several hundred stations. Thus, 
the total number of variables computed approaches nearly 
one-million for each case. In this section results for two 
typical cases will be presented. These results consist 
of hydrogen concentration, temperature and pressure 
distributions at three axial locations and at two locations 
transverse to the main flow direction. 

9 . 2- 2 Definition of cases 


One case for normal injection (Case 1) and two cases for 
parallel injection (Cases 9 and 7) are described in this 
chapter. The geometrical configuration of these cases 
are shown in Figures l, 2 and 3 respectively. Table 6 
defines the flow and thermal properties. 


79 



TABLE 6;. DEFINITION OF PROPERTIES 


Location 


Main Stream 
conditions 


Jet 

conditions 


Boundary 


Property (units) 


Flow speed (m/s) 

Temperature (°K) 

Mass fraction Ng 

Mass fraction Og 

Mass fraction HgO 
2 

Pressure (N/m ) 


Flow speed (m/s) 

Temperature (*^K) 

Mass fraction Hg 
2 

Pressure N/m 
Mass flow Kg/sec 



Case 1 


675 

70.6 

.7676 

.2325 

0.0 

8720.0 



1585.2 

1178.6 

.487 

.263 

.25 

179300.0 


Case 7 


1585.2 

1178.6 

0.487 

0.263 

0.250 

179300.0 


1210.0 

2039.4 

2039.4 . 

250.0 

150.3 

150.3 

1.0 

1.0 

1.0 

212000.0 

179300.0 

179300.0 

.000165 

.0109 

0.0109 

J 

Ducted* 

Free 

Ducted* 
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9.3 Graphs of Results 


All data are presented for planes normal to the wall and 
parallel with the main flow direction. One of these planes 
is located in the transverse direction directly at the ' 
hydrogen jet centerline and the other is located between 
jets (see Figures 1 and 2). The comparison of profiles at 
a given plane at various main flow stations Indicates the 
rate of movement of hydrogen in a direction normal to the 
flow. * 

V 

Figures 13 throughl5 are for normal injection (Case 1). 
Figure 13 shows the distribution of hydrogen in any form at 
three axial locations (.118m, .216m and .246m). The first 
profile is before the injection point (.186m) and therefore 
no hydrogen is present. At the .216m location, the 
hydrogen concentration levels (less than .01) indicate that 
substantial mixing is primarily in the y direction. (The 
close spacing of the adjacent injectors limits x direction 
mixing). Examination of the results at the .246m station 
shows continuing movement of hydrogen upward with an 
attendant smoothing of the profile, as would be expected. 

The temperature profiles of Figure 14 indicate that after 
injection and combustion that little transverse thermal 
gradient exists. This is partly due to the fact that the 
flame front is at the outer layer of the hydrogen zone 
which tends to diffuse the temperature in the transverse 
direction. Also it is seen that the reaction of hydrogen 
moves the thermal boundary layer outward from the plate 
with movement in the main flow direction. 

The pressure distribution for Case 1 is shown in Figure 15 
It exhibits the qualitative behaviour one would expect. 
There is little transverse pressure difference due to low 
velocities in that direction. A pressure spike can be 
seen downstream of the jet indicating a shock. The 
pressure below the shock is essentially uniform at a 
higher value than the free stream. It must be emphasized 
that the pressure distribution in the jet region is subject 
to error due to the fact that the flow was forced to be 
parabolic when in actuality there is some upstream effect 
caused by the jet. 

Figures 16 through 18 are results for Case 9 which is 
parallel injection. The diameter of the jets and jet 
spacing are much larger than for Case 1. The profiles at 
injection (z=0) show that there is little mixing of 
hydrogen between jets. However, as the flow moves forward 
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FIGURE 18: PRESSURE DISTRIBUTION FOR CASE 9 



the mixing increases such that the concentration of hydrogen 
between jets is roughly 50% of the hydrogen concentration in 
line with the jet at zp.36 meter. However, these figures 
show that in the normal direction the movement of hydrogen 
upward between jets is not very effective. This is due to 
the low normal velocities for the parallel injection case. 

The temperature distribution in Figure 17 indicates a higher 
temperature between jets even though the hydrogen 
concentration is lower. This is because the combustion 
occurs at the fringe of the hydrogen zone. The combustion 
zone upper surface is indicated by the "spikes" in the 
temperature profile. This combustion is seen to grow with 
movement downstream. The pressure distribution shown in 
Figure 18 does not indicate any shocks. At the injection 
station ( 2 = 0 ), the flow expands in the low density region. 
However, at downstream locations the pressure becomes 
uniform at approximately the free stream value. 

9.3-1 Parallel Injection into a duct with expanding top wall 

This is the case chosen to demonstrate the validity of the 
SHIP program. 

Profiles of 19 variables described earlier were computed for 
a grid of 12 grid lines transverse to the flow and 12 grid 
lines normal to the flow for a total of 144 points at each 
station in the main flow direction. Computations were made 
up to 2=0. 50m (see Figure 3) for 150 stations. In this 
section the results thus achieved will be presented. They 
consist again of hydrogen concentration, temperature and 
pressure distribution at four axial locations and at two 
locations transverse to the main flow direction. (Figures 19, 
20 and 21,) . 

The concentration profiles at injection (z=0) (Figure 19) 
show again that there is little mixing of hydrogen between 
jets. Downstream the same arguments as for Case 9 above 
still hold true. The temperature distribution in Figure 20. 
indicates again a higher temperature between jets even though 
the hydrogen concentration is lower. The pressure distribution 
(Figure 21) does not indicate any shocks. At the injection 
station (z=0), the flow expands in the low density region. 

At downstream locations the pressure becomes uniform at 
approximately the free stream value up to the point where 
expansion of the domain starts. Downstream that point the 
pressure is again uniform, but at lower values as we proceed 
further downstream. 
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CONCLUSIONS AND RECOMMENDATIONS 


Qualitatively the results obtained by the method developed 
herein appear to be correct. Further, comparison with cold 
flow data also show excellent agreement. The computational 
time for these results is also very good. 

The assumption of parabolic flow in the injection region 
needs improvement because both data and calculations made 
in this work indicate that the pressure distribution is 
affected upstream by the jet. Recirculation, however, is 
probably not important. What is needed is a technique 
to allow the pressure to be calculated via an elliptic 
procedure and the velocity by a parabolic one. Such a 
technique is advantageous over fully elliptic procedure 
because storage and computation time are much lower. 

The SHIP code has considerable potential for further 
development and exploitation. The suggestions made here 
are by no means exhaustive, but are merely intended to 
be pointers. Two distinct (but related) possibilities 
exist . 

Firstly, the code can be extended to handle other, more 
complicated situations by suitable modifications. 

Secondly, it can be varied to produce new codes capable 
of solving new problems. At the present time the 
following suggestions seem worth recording. More 
improvements may be introduced into the SHIP code to increase 
its efficiency. As a matter of fact a number of latest 
advances are now available and may be transferred to the SHIP 
code. They include major reprogramming for improving the 
numerical algorithm (e.g. use of a different differencing 
scheme) and making the code more readily understandable. 

A disadvantage of the present version of the code is that 
both velocities and pressures are stored at the same axial 
position (i.e. they are not staggered in the predominant 
flow direction). Staggering the velocity locations even 
in that direction and using different density interpolation 
formulae, will have several advantages; apart from others 
this practice has also the distinct merit that it blends 
smoothly with the elliptic and partially-parabolic ones. 
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APPENDIX A 


The Chemical Equilibrium Model 

The main features of the equilibrium chemistry model used 
to predict the properties in a hydrogen - air flame are 
described below. 


Four equilibrium reactions are assimed as follows: 


H + 

H t 

»2 

(A-1) 

0 + 

0 t 

^2 

(A-2) 

H + 

^3 
OH t 


(A- 3)' 

0 + 

H t 

OH 

(A-4) 

The 

six species 

involved in these reactions are considered 


to be present with nitrogen which is inert. In developing 
the equations to predict the equilibrium concentration of 
the species, two quantities are defined, namely: 


II 

CO 

+ m^ 

. *0 


(A-5) 


^2° ^OH 

f = m„ 
^2 

H- mjj 


H 2 O Wqjj '"oh 

(A-6) 


where F is the total mass fraction of oxygen in any form 
and f is the total mass fraction of hydrogen in any form. 
Since the molecular weight of the various oxygen species 
is approximately equal to that of nitrogen, it is assumed 
that the coefficient of turbulent diffusion of the chemical 
species are equal to each other at every point in the flow. 
A well known consequence is that f is linearly related to 
nijj and the constants in the relation can be determined 

from the boundary conditions: 



1 - f 


(A-7) 


2 , air 


2 . air 


The total mass fraction of all elemental species must be 
unity which gives, 


^ “Ng 1 

Also, 


(A-8) 


F = q (1 - f) 
with q = m« = 0.232 

^2 , ai r 


(A-9) 


Therefore, if f is known, F can be determined by equation (A-9) 


From thermodynamic considerations the equilibrium constant 

K , for the reaction: 

P 

K 

P 

aA + bB cC is defined by. 


c 

K 

a b 
A 


c-a-b 


(A-10) 


where x stands for concentration and the pressure p is in 
atmospheres. For each of the four reactions (A-1 ) to (A-4) in 
the present model c-a-b = -1. Expressing the concentrations 
in terms of mass fractions by noting that. 


Wi 

*"i = ¥“ ^i 


(A-11) 


where W is the molecular weight of the mixture we get: 


”A B A "^B 


(A-12) 


Thus the equilibrium equations for the reactions (A-1) to 
(A-4) can be written: 


(A-13) 
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(A-14) 


I 



m2 

0 

"’h’^oh 

"’oh 


(A-15) 


(A-16) 


The condition of equilibrium is expressed using four 
equilibrium constants for the four chemical reactions, 
thermodynamic equilibrium prevails, the K’s take values 
which depend upon temperature alone. In the above set of 


If 


equations there are seven unknowns (m 


H, 


m. 


m 




m 


m. 


0’ "OH 


m. 


HgO and F; f being given) and seven equations (A-5, A-6, A-8,- 


A-13, A-14, A-15, A-16), The ^problem is therefore soluble. 


The remaining discussion defines the solution procedure. 
Derivation of Solution Procedure 


The present procedure is based on the reduction of the 
number of variables under consideration. Two equations 
are derived as follows: 

Equations (A-13 - A-16) are solved to give the relationships: 


m. 


'H 










t/tcC 


(A-17) 


(A-18) 
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(A-19) 


%0 


■% h 


k[ 




“‘oh 





1 ^2 




(A-20) 


It is convenient to define the following parameters*: 





B 



C 


D 


k; k; 


! 




f ? 


K. 



I 



(A-21) 


V (A-22) 


(A-23) 


(A-24) 


♦These depend upon temperature alone » therefore » they can 
be tabulated right at the start. 
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Using equations (A-17 to A-20) to eliminate mQ, q, and 

2 

equations (A-6) and (A-6) gives: 

f = mjj (1+ic }/in^ ) + ^A + — D /m^ ')(A-26) 

2 2 ' 2 2 

Where the definitions given by equations (A-21 to A-24) have been 
used. 

These equations have the form of a quadratic equation: 

— 2 - - 
au + bu + c = 0 

with solution 


u 


-2 c 



4ac 


(A-27) 


Note that a and b are always positive and c is always 
negative. Since u > 0, the physically meaningful root is 
the one with the positive sign. This particular form of 
quadratic expression is chosen since it does not require 
subtraction and gives greater precision. Using equation 
(A-27) to express the solution of equations (A-25) and (A-26) 
gives : 
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f 




A * 17 ® -S, 


7(1+ i )f (A-29) 


Equations (A-28)and(A-29)contain the unknowns and in„ . 

''^2 "2 

Before detailing the solution procedure several properties of 
these equations will be discussed. Table A-1 gives the values 
of k| , Ko, Kg, , A, B, C and D for the temperatures from 
100 to OOOOK'^at a constant product of 2.5 (corresponding 
approximately to 1/10 atmosphere). Note that while the 
individual constants vary from lesS' than 1 to 10^*^®, the 
groups appearing in equations(A-28) and( A-29 ) vary much less 
and can be easily processed on a digital computer. Second 
it is to be observed that the equations have been derived 
so that a 'physical' relationship has been established 


H, 


0 , 


and f , i . e . as F 0 , m^ 
should also be observed ?hat m. 


between m„ and F, and m. 

^2 

and as f ->• 0 , mjj ->• 0 • 

and m„ can neve^ be less than zero or greater than 1. 

“2 

Finally, neither equation becomes indefinite as m„ or m^ 

approach zero. Thus, the equations. are well-behaved and 
can be readily solved for wide variations in temperature 
and pressure. The solution procedure is described next. 


“0, 


value of m„ is guessed in the following wav, 

.2 


If the value 


at the equivalent upstream station is known, then that value 
is used; otherwise, a value of zero will always lead to 
a converged solution. This assumed value of m„ is substituted 

into equation (A-28)which yields mQ . Then the computed 

value of m-, is substituted into equationC A-29)which allows 
^■'2 

calculation of a new value of m 


H, 


The assumed and calculated 


values of m„ .are compared. If ?hese values differ by more 


^2 

than a specified convergence criterion, the calculated value 
of m„ is taken as the assumed value and the process described 
”2 

above is repeated until convergence is obtained. The 
behaviour of this solution technique is shown in Figure A-1. 
This figure shows with the broken line a plot of mj^ 
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FIGURE A-1: GRAPHICAL PLOT OF TRI AL-AND-ERROR SOLUTION 


calculated versus m„ assumed. The correct value is achieved 

when the two values are equal . The locus of points for this 
situation is a straight line with a slope of unity. If 
nijj assumed is less than the correct value, the figure shows 

that the calculated value will always be larger. Thus, 
when this calculated value of m„ is taken as the assumed 

^2 

value, the resulting newly calculated value will be closer 
to the correct one. The same argument can be made to 
show that if the initial assumed value of m„ , is too large, 

Hg 

the iteration process will again cause convergence to the 
correct one. Great precision can be obtained with this 
method. The convergence criterion for the calculation is: 

jmjj calculated - nij^ assumed | 0.001 m^ calculated 

2 2 2 

Greater precision can be easily achieved by a stricter 
convergence criterion. However, the present one is 
sufficient for most practical calculations. 

It may be concluded that the method offers a reliable, 
simple and 'extremely fast technique for solving the 
equilibrium equations arising from the chemical model 
considered in this work. 
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TABLE A-1. NUMERICAL VALUES OF EQUILIBRIUM CONSTANTS 


Temp 

k: 

1 

^2 


^4 

200 

5.0+100 

3.1E+99 

2.6+100 

2.6+100 

400 

2.8E+52 

2.7E+58 

3.2E+51 

6.2E+59 

600 

2.3E+33 

4.3E+ 36 

5.8E+32 

9.0E+37 

800 

5.9E+23 

4.9E+25 

2.2E+23 

9 . 3E+26 

1000 

9.7E+17 

1.2E+19 

4.5E+17 

2.2E+20 

1200 

1.2E+14 

5.0E+14 

7.0E+13 

8.2E+15 

1400 

2.1E+11 

3.5E+11 

1.2E+11 

5.4E+12 

1600 

1. 7E+09 

l.lE+09 

l.lE+09 

2.2E+10 

1800 

3.9E+07 

2.1E+07 

2.7E+07 

3.0E+08 

2000 

1.9E+06 

7.0E+05 

1.4E+06 

9.8E+06 

2200 

1.5E+05 

4.3E+04 

1.2E+05 

5.8E+05 

2400 

1.9E+04 ■ 

4.2E+03 

1.6E+04 

5-5E+04 

2600 

3.4E+03 

5.8E+02 

2.8E+03 

7.6E+03 

2800 

7.5E+02 

l.OE+02 

6.4E+02 

1.3E+03 

3000 

2.0E+02 

2.4E+01 

1.74+02 

3.1E+02 

3200 

6.3E+01 

6.8E+00 

5.7E+01 

8.5E+01 

3400 

2.3E+01 

2.1E+00 

2.1E+01 

2.7E+01 

3600 

9. 3E+01 

7.9E-01 

8.6E+00 

9.8E+00 

3800 

4.1E+00 

3.2E-01 

3-8E+00 

3.9E+00 

4000 

1.9E+00 

1.4E-01 

1.8E+00 

1.7E+00 

4200 

l.OE+00 

6.8E-02 

9.8E-01 

8.3E-01 

4400 

5.5E-01 

3.5E-02 

5.4E-01 

4.2E-01 

4600 

3.2E-01 

1.8E-02 

3.1E-01 

2.2E-01 

4800 

1.9E-01 

l.OE-02 

1.9E-01 

1.3E-01 

5000 

1.2E-01 

6.4E-03 

1.2E-01 

7.7E-02 

5200 

7.9E-02 

4o0E-03 

7o8E-02 

4.8E-02 

5400 

5.3E-02 

2»5E-03 

5.3E-02 

3.1E-02 

5600 

3.6E-02 

1.7E-03 

3.6E-02 

2«0E-02 

5800 

2.6E-02 

iaE-03 

2o6E-02 

1.4E-02 

6000 

1.8E-02 

8.1E-04 

1.8E-02 

9.8E-03 


A 

B 

C 

5 

4.4E-51 

1.7E-50 

2.42E+50 

1 .95E+00 

5.9E-27 

6.0E-30 

4 .30E+29 

2.19E+04 

2.0E-17 

4.7E-19 

1.09E+19 

8.46E+02 

1.2E-12 

1.4E-13 

4.95E+13 

1.56E+02 

l.OE-09 

2.7E-10 

2.94E+10 

5.94E+01 

8. 7E-08 

4.4E-08 

2.14E+08 

3.14E+01 

2:iE-06 

1.6E-06 

5.22E+06 

1.81E+01 

2.4E-05 

2.5E-05 

3.67E+05 

1. 32E+01 

1.5E-04 

2.1E-04 

4.53E+04 

9.45E+00 

7.2E-04 

l.lE-03 

8„63E+03 

7.76E+00 

2.5E-03 

4.8E-03 

2.24E+03 

6.96E+00 

7.0E-03 

1.5E-02 

0.71E+03 

5. 77E+00 

1.7E-02 

4aE-02 

2. 60E+12 

5 . 30E+00 

3.6E-02 

9.6E-02 

l.llE+02 

4.49E+00 

7-0E-02 

2.0E-01 

5.37E+01 

4.34E+00 

1.2E-01 

3.8E-01 

2.95E+01 

3.87E+00 

2.0E-01 

6.7E-01 

1.70E+01 

3.62E+00 

3.2E-01 

l.lE+00 

1.02E+01 

3.45E+00 

4.9E-01 

1.7E+00 

6.38E+00 

3.25E+00 

7.0E-01 

2.6E+00 

4. 30E+00 

3.09E+00 

9.8E-01 

3.8E+00 

3.12E+00 

3.09E+00 

1.3E+00 

5.3E+00 

2.20E+00 

2.89E+00 

1.7E+00 

7.2E+00 

1.59E+00 

2.69E+00 

2.2E+00 

9.6E+00 

1.30E+00 

2. 74E+00 

2.8E+00 

1.2E+01 

9.62E-01 

2.59E+00 

3.5E+01 

1.5E+01 

7.49E-01 

2.52E+00 

4.3E+00 

1.9E+01 

S.20E-01 

2.53E+00 

5.2E+00 

2o4E+01 

4.85E-01 

2.49E+00 

6.1E+00 

2o9E+01 

4.22E-01 

2.47E+00 

7.2E+00 

3.5E+01 

3.44E-01 

2.47E+00 
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APPENDIX B 


Free-Stream Boundary Conditions in 
Supersonic Plows 

(a) Consider the effect of a pressure wave of angle P on 
a stream of velocity Woo (supposing that v » 0 and 
« 0 ). 



The v-momentum equation for the domain ABCD (where 
AB and CD are stream surfaces) is; 

= (p„ -p) / tans. (B-1) 

tan S = (Wg - w„)/(v 2 - v^) 

■- The wave angle can be related to the Mach number by: 

tans « 1// wf - 1 ( 3 _ 2 ) 

so that: 

p - Poo r~2 ~ 

V = s/vlI . 1 (B-3) 

P w 

OO CO 


(b) The energy equation for a perfect gas gives: 

. uj , Y P-K vP 

2 (Y-l)P^j “2 (V-1) P^ 


(B-4) 
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where Uj, denotes the resultant velocity at the 
boundary. For an isehtropic flow. 



whence: 



CO 


, (B-5) 
, (B-6) 


SO that: 


(Pb/^co> 




(1 - 1/T) p. 


Substitution of (B-7) into (B-4) gives: 



and of course 



(B-7) 


(B-8) 


(B-9) 
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APPENDIX C 


List- of FORTRAN Variables 


Notes 

1. ■ Star rating 

This list classifies variables as having general 
importance or being merely of local significance. 

The latter variety is marked with a dash (-). In 
the former variety, each variable is given a star 
rating of one (*), two (**) or three (***) stars. 

This indicates that the variable is appropriate to 
the problem-dependent (*), physical -modelling (**), 
or the main-machinery (***) part of the program. 

2. Subscripted variables 

Subscripts to variables are shown only • when necessary, 
in order to explain the meaning of the variable. In 
general therefore, the subscript status of each 
variable should be understood from the column headed 
NATURE . 

3 . Symbols 

Wherever possible, the symbols corresponding to 
FORTRAN variables, which are used in the text, are 
also listed. 
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VARIABLE 

NATURE 

1 

A(I) 

REAL, ARRAY 

2 

ADIN 

REAL, 

3 

ADS(IJ) 

REAL, ARRAY 

4 

AGEOM 

REAL, 

5 

AK 

REAL, 

6 

AKE(IJ) 

REAL, ARRAY 

■ 

AKFAC 

REAL, 


AKIN 

REAL, 

1 

ALX,ALXM, 
ALXP , ALY 
ALYM , ALYP 
ALZ,ALZ1 

REAL, 

10 

AMUREE 

REAL, 

11 

AMUT(IJ) 

REAL, ARRAY 

12 

ARAT 

REAL, 

13 

AREA 

REAL, 

14 

AXM,AXP, 
! AYM,AYP, 
AZ 

I 

REAL, ARRAY 

15 

B 

■ 

REAL, ARRAY 

16 

BXE 

REAL, 

m 

BXEU 

REAL, 

18 

BYN 

REAL, 


BYNU 

REAL, 


CD* 

REAL, 


STAR 

RATING 

SYMBOL 

MEANING 



Transformed coefficient 
in TDMA operation 



Dissipation rate at inlet 


e 

Dissipation rate 

* 

- 

Factor in grid expansion 


K 

Mixing-length constant 

*** 

k 

Turbulence energy 

— 


Factor relating k to 
mean motion energy 

* 


Inlet kinetic energy 


C 

^e, 

^2 etc 

Quantities representing 
the convective 
coefficients in the 
finite-difference 
equations 

♦ 


Reference laminar 
viscosity 



Turbulent viscosity 

* 


Ratio of upstream to 
downstream areas 

- 


Area of control-volume 
face 


^W’^E' 

^S’^N, 

^P,U 

Coefficient of finite- 
difference equations 



Transformed 
coefficient in TDMA 
operation 


"^E 

Width of integration 
plane in the ^-direction 

* 


Upstream value of BXE 

** 

Arij^ 

Width of integration 
plane in the n- 
direction 

★ 


Upstream value of BYN 


Cd 

Constant in turbulence 
model 











NO 

VARIABLE 

NATURE 

STAR 

RATING 

SYMBOL 

. 

MEANING 

21 

CDQR,CDBT, 

CDTQ 

REAL, 

’ 

** 


Roots of CD 

22 

CP 

REAL, 

* 

Cp. 

Specific heat 

23 

CPDCV 

REAL 

* 


Specific heat at 
constant pressure 
divided by specific 
heat at constant 
volume 

24 

CX, CXP, 

CY, CYP, 

CZ , CZP 

REAL, CY only 
as ARRAY 

■ 


Temporary storage for 

finite-difference 

coefficients 

25 

C1,C2 

REAL, 

* 


Constants in 
'turbulence model 

26 

CFS etc 

REAL, 

♦ 


Local shear stress 
at South wall 

27 

CFSAV,etc 

REAL, 

* 


Local normalised 
shear stress at South 
wall 

28 

DBXEDZ , 
DBYNDZ 

REAL, 


A(i5E>. 

dz 

dz 

Rates of growths of 
BXE, BYN with respect 
to dz 

1 

,29 

DEN • 

REAL, 


P 

Reference density of 
fluid 

30 

DHED 

REAL, 

* 


Dynamic Head 

31 

DHEDIN 

REAL, 

♦ 


Value of DHED at inlet 

32 

DPDZ . 

REAL, 


3z 

Pressure-gradient in 
the c-direction 

33 

DPDZU 

REAL, 

* 


Upstream value of 
DPDZ 

34 

DU,DV,DW 

REAL, ARRAY 

*** 

Pressure-difference 
coefficients for the 
U,V & W velocities 

35 

DXWDZ, 

DYSDZ 

REAL 

** 


Slopes of W & S 
boundaries 

r 

DZ ' 

REAL, 


AS 

Forward step-size 
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■' '"I 

NO 

VARIABLE 

NATURE 

37 

DZU 

REAL, 

38 

EE 

REAL, 

39 

ETA 

REAL, ARRAY 

40 

EX 

REAL, 

41 

F 

REAL, ARRAY 

42 

FLOWIN 

REAL, 

43 

FLOINJ 
FLO I NT 

REAL, ARRAY 

44 

FLUXE , 
FLUXW, 
FLUXN, 
FLUXS 

REAL, ARRAY 

45 

FRA, 

FRAM 

REAL, 

46 

PXM,FXP, 

FYM,FYP 

REAL, ARRAY 

47 

GAM 

REAL, ARRAY 

48 

GAMA 

REAL, 

49 

GAME, 

GAMW, 

GAMS, 

GAMN 

REAL, ARRAY 


STAR SYMBOL 
RATING 


MEANING 


Upstream value of DZ 

Constant in the law- 
of-the-wall 

Non-dimensional y 
coordinate 

Factor by which 
forward step size is 
incremented 

Store for dependent- 
variable values 

Rate of mass inflow 
into the calculation 
domain 

Mass flow rate 
injected by a jet 
through the south and 
north wall respectively 

Fluxes of dependent 
variables on E,W,N,S 
boundaries 


Fraction of boundary- 
layer width used in 
calculating forward 
step size, fram is the 
maximum value of FRA 

Interpolation factors 
to indicate distance 
of a pressure node 
from the neighbouring 
velocity locations 

r. Store for diffusion 

coefficients 

Local value of GAM at 
locations where it is 
not stoi’ed 

r. Boundary values of GAM 

^ ’ to compute corresponding 

FLUX'S 
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NO 

VARIABLE 

NATURE 

STAR 

RATING 

SYMBOL 

^ MEANING '* 

50 

GASCON 

REAL, 

* 

R 

Universal gas constant 

51 

• GREAT 

REAL, 



Large number, used to 
fix values to desired 
levels and in limiting 
overflows; 10^® 

52 

GX,GY 

GZ 

REAL, ARRAY 

** 


Mass velocities in the 
n and C-directions 

53 

HO 

REAL, ARRAY 

♦ 


Enthalpy of formation 

' 54 

I 

INTEGER, 


Index indicating 
position in the 
C-direction 

55 

I C JUMP 

INTEGER , 

* 


I n dex , controlling 
calculation of wall 
fluxes for purposes 
of printout only 

56 

IINJ, 

IINJT 

INTEGER, 

* 


Indices to control 
entry to INJMOD and 
INJMOT, respectively 

57 

IJ,IJE 
IJN,IJNW, 
US, USE, 
IJW 




Computed subscripts 
to replace ( I , J ) , ( I+l , J) 
(I,J+1),(I-1,J+1),(I-1, ' 

J),(I+1,J-1), (I-1,J) 

58 

IJREF 

INTEGER, 

*** 


Location of grid node 
that is used as the 
reference for relative 
pressure 

59 

IJT 

INTEGER, 


■ 

Calculable subscript 
for T (*IJ+NFM(NVT)) 

60 

I MAX 

INTEGER, 

* + ♦ 


Maximum value of I , ; 
for which storage j 

locations are provided,, 
in the program ' 

61 

INJSTR, 

INJTOP 

INTEGER, 

* 


Value of ISTEP at^ , 
j njection : location , , 
for South and North; 
wall, resjbectively 

62 

IPRINT 

(NV) 

INTEGER, ARRAY 



Values of a variable NV, 
printed only if the 
corresponding IPRINT 
is unity 
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NATURE 


STAR SYMBOL 
RATING 


MEANING 


I SOLVE 
(NV) 


INTEGER, ARRAY 


64 I STEP INTEGER, 

65 ISTR INTEGER, 


66 ISWP 


67 IXY 


68 J 


70 JMAX 


JSMAX 


JSTR 


JSW? 


INTEGER, 


INTEGER, 


INTEGER, 


69 JM(J) INTEGER, ARRAY 


INTEGER, 


INTEGER, 

INTEGER, 


INTEGER, 



Index denoting whether 
the finite-difference 
equations of variable 
NV are solved 
(I SOLVE (NV) = 1) or not 
(ISOLVE(NV) = 0) 

Counter for foiward 
steps 

Value of I for the 
first internal storage 
location for a given 
variable 

Index denoting the 
direction of sweep 
while performing the 
TDMA traverse in the 
^-direction 

Integer denoting 
whether the 
C-direction or the 
n- direct ion TDMA 
traverse is performed 
first 

Index denoting the 
position in the 
n-direction 

(J-1)*IMAX ; used for 
subscript calculation 

Maximum value of J, 
for which storage 
locations are 
provided 

J-locatlon of SMAX 

Value of J denoting 
the first internal i 


location for a given 
variable 

Index denoting the 
direction of sweep 
while performing the 
TDMA traverse in the j 
n-direction 








NO 


VARIABLE 


NATURE 

• STAR 
RATINC 

SYMBOL 

t 

MEANING 

INTEGER, 



Indices, denoting the 
nature of the E.W.N 
and S boundaries; 

= 1 for wall 
boundaries; = 2 for 
symmetry planes; = 3 
for 'free* boundaries 

INTEGER, 



Number of grid lines 
minus one in the 
f;-direction ; 

INTEGER, 

* 


The LAST STEP i.e. 
the maximum number 
that ISTEP can attain 

INTEGER, 



The number of main 
control volumes in the 
^-direction 

INTEGER, 

3k 9k 9k 


Number of grid points 
in the ^-direction 

INTEGER, 

9k 9k 9k 

I 

Number of grid lines 
minus one in the 
n-direction 

INTEGER, 

9k 3k 9k 


Number of control 
volumes in the 
n-direction 

INTEGER, 

9k 9k 9k 


Number of grid points 
in the n-direction 

INTERGER, ARRAY 

9k 9k 9k 


(NV-1)*IMAX*JMAX; used 
for subscript 
calculation 

INTEGER, 

9k 


Maximum number of 
variables F, for which 
arrangement for 
Printout is available 

INTEGER, 

9k 9k 


‘ Number of iterations 
on variables <j> (i.e. 
F) other than 
velocities 

INTEGER, 

9k 9k 


Number of iterations 
on velocity components 

INTEGER, 

9k 9k 9k 


Index identifying GAM, 
the effective diffusion 
coefficient 
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75 

76 

77 

78 
80 

81 

82 

83 

84 

85 

86 

87 


KBCE, 

KBCW, 

•KBCN, 

KBCS 


LASTEP 

LCV 

LPl 

M 

MOV 

MPl 

NFM(NV) 

NFPMAX 

NITERF 

NITERM 

NMIJ 
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NO 

VARIABLE 

NATURE 

STAR 

RATINC 

SYMBOL 

r 

88 

NNV 

INTEGER, 



89 

NH 

INTEGER, 



90 

■ NH2 

INTEGER, 



91 

NH20 

INTEGER, 



92 

NN2 

INTEGER, 



93 

NO 

INTEGER, 



94 

NOH 

INTEGER, 



95 

N02 

INTEGER, 



96 

NPJUMP 

INTEGER , 

* 


97 

! 

NPP 

INTEGER, 

*** 


98 

NRHO 

INTEGER, 

*** 


99 

NSWP 

(NV) 

INTEGER, ARRA^ 

*** 


100 

NV 

INTEGER, 




MEANING 


Number of variables 
for which storage in 
the F array is provided 

Identifier of atomic 
hydrogen concentration 

Identifier of molecular 
hydrogen concentration 

Identifier of water 
vapor concentration 

Identifier of molecular 
nitrogen concentration 

Identifier of atomic 
oxygen concentration 

Identifier of OH 
concentration 

Identifier of 
molecular oxygen 
concentration 

Intervals of ISTEP 
after which a 
printout of main 
variables is obtained 

Index identifying p' 
in the F-array 

Index identifying p 
(density) 

Number of pairs of 
TDMA sweeps for 
variable (NV) 

Serial number of any 
variable <|) stored in 
the F-array 
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NATURE 


STAR SYMBOL 
RATING 


101 • N VF , N VD INTE GER , 

NVK,NVH 


102 NVP,NVT, INTEGER 
NVU.NVV, 

NVW 



104 PI 

105 PIN 

106 PJAY 

107 pp 

108 PR(NV) 

109 PREF 

110 PRLAM 
(NV) 

111 RELAX 
(NV) 

112 


REAL, ARRAY 
REAL, 

REAL, 

REAL, 

REAL, ARRAY 
REAL, ARRAY 

REAL, 

REAL, ARRAY 
REAL, ARRAY 



RE TRAN REALi 


MEANING 


Identifiers for hydrogen 
concentration, ( in any 
form) , dissipation rate 
e, turbulence energy 
k, and enthalpy 
n, in F 

Identifiers for pressure 
p, temperature T, and 
velocity components u, 

V and w, respectively 

Static pressure p 

Constant 

Inlet value of p 

Resistance of laminar 
sublayer 

Pressure correction 

Effective Prandtl/ 
Schmidt number 

Pressure at reference 
grid node 

Laminar Prandtl/Schmidt 
number 

Relaxation factor for 
NV 

Transition 
Reynolds number 


113 RHO(IJ) REAL, ARRAY 

114 ROAR REAL, 

115 SMALL REAL, 


Density 

Density and control . 
volume face area 
product 

Sitial 1 number used to 
trap the possibility 
of divisions by zero ; 
10-30 used 
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NO 

VARIABLE 

NATURE 

STAR 

RATINi 

SYMB 

G 

116 

SP 

REAL, ARRAY 


Sp 

117 

STANS 

etc 

REAL 

♦ 


118 

STANSM 

etc 

REAL 

* 


119 

SU 

REAL, ARRAY 

*** 

Su 

120 

TIN 

REAL, 

* 


121 

TITLE 

(NV) 

REAL, ARRAY 



122 

TINJ, 

TINJT 

REAL 

♦ 


123 

TWAL 

REAL, 

* 


124 

TX,TY 

REAL, 

- 

T?t' 

125 

U(IJ) 

REAL, ARRAY 

** * 

u 

126 

1 

TJIN 

REAL, 



127 

V(IJ) 

REAL , ARRAY 


V 

128 

VIN 

REAL, 

♦ 


129 

VI NJ, 
VINJT 

REAL, 

* 


130 

VOL 

1 

REAL, 

- 



MEANING 


One part of the 
linearised source 
term 

Local Stanton numbers 
■ at South wall, etc. 

Mean Stanton 

number at South wall, etc 

A part of the 
linearised source 
term 

Inlet value of 
Temperature T 

Array storing 
36-character alpha- 
numeric names of 
members of F 
Inlet jet temperature 
at South and North 
walls respectively 

Temperature of wall 
boundaries- 

Temporary store for 
the diffusive coefficients 
Velocity component in 
the 5-direction 

Inlet value of u 

Velocity component 
in the n-direction 

Inlet value of v 

Velocity of jet 
injection at South 
and North walls, 
respectively 

Volume of control 
volume, surx’ounding 
a grid node 








NO I VARIABLE NATURE 


131 W(IJ) 


132 WEAR 


STAR SYMBOL 
RATING 


133 WIN 

134 WM 

135 ZD 


136 ZLAST 


137 ZU 

138 ZETA 

139 ZINJ, 
ZINJT 


140 ZRE 




MEANING 


Velocity components 
in the ^-direction 

Mean value of w across 
each plane of 
calculation 

Inlet value of w 

Molecular weight 

Downstream value of C 
at which calculations 
are currently being 
performed 

Last value of c a-t 
which calculations are 
to be performed 

Upstream value of C 

Non-dimensional x 
coordinate 

^-location of 
injection at South 
and North walls, 
respectively 

Values of z for which 
printout is desired 
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APPENDIX D 


Listing of SHIP . 

The subroutines and their entry points are listed in the 
following order. 


1. 

BLOCK 

2. 

MAIN 

3. 

ALLMOD 


• 

BEGIN 


• 

GAMOD 

- 

• 

GEOMOD 


O 

UPSTRM 


• 

SOMOD 



INJMOD 


• 

INJMOT 

4. 

AUX 


• 

DENSTY 


• 

GAMMA 


• 

SOURCE 


• 

VISCOS 


• 

SPECIE 

5. 

PRINT 

6. 

SOLVE 

7. 

STRIDA 


• 

STRIDO 


• 

STRIDl 


• 

STRID2 

8. 

STRIDE 


• 

STRID3 


• 

STRID4 
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